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ABSTRACT 


The  calculations  of  heat  losses  and  temperature  fields  for  directly  buried 
conduit  heat  distribution  systems  were  performed  using  finite  element 
method.  The  finite  element  analysis  solved  two-dimensional,  steady-state 
heat  transfer  problems  involving  two  insulated  parallel  pipes  which  were 
encased  in  the  same  conduit  casing  and  in  separate  casings,  and  the 
surrounding  earth.  Descriptions  of  the  theoretical  basis,  computational 
scheme,  and  the  data  input  and  outputs  of  the  developed  computer  programs 
are  presented.  Numerical  calculations  were  carried  out  for  predicting  the 
temperature  distributions  within  the  existing  high  temperature  hot  water 
distribution  system  with  two  insulated  pipes  covered  in  the  same  metallic 
conduit  and  the  surrounding  soil.  The  predicted  results  generally  agree 
with  the  experimental  data  obtained  at  the  test  site.  The  deviations 
between  the  predicted  and  measured  values  are  found  to  range  from  0 to  17 
percent  with  an  average  of  6 percent.  The  rates  of  heat  loss  from  two 
insulated  pipes  encased  in  separate  conduits  were  calculated  for  different 
pipe  sizes,  fluid  temperatures  and  insulation  thicknesses.  The  results 
were  compared  with  the  predictions  from  a steady-state,  one -dimensional 
radial  heat  conduction  model.  The  discrepancies  between  finite  element  and 
radial  conduction  models  in  pipe  heat  loss  values  are  discussed. 

Key  words:  Computer  program,  direct  burial,  district  heating  and  cooling, 
finite  element  method,  heat  loss,  heat  transfer,  underground  heat  distribution 
system. 
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1 . Introduction 


District  heating  and  cooling  systems  can  deliver  hot  water  and  chilled 
water  economically  and  efficiently,  from  the  central  plant  to  the  end  users 
through  insulated  piping  networks.  A centralized  system  is  often  more 
economical  to  produce  heating  and  cooling  energy  since  it  cost  less  to  use 
local  energy  sources  or  to  burn  low-cost  fuels  such  as  coal  and  municipal 
solid  waste.  In  addition,  a more  efficient  delivery  of  energy  can  be 
obtained  with  district  heating  and  cooling  systems  due  to  reduced  operating 
and  maintenance  costs  and  improved  pipeline  design  and  construction  methods 
compared  to  individual  local  building  plants. 

Military  institutions  and  facilities  maintain  approximately  6000  miles  of 
heat  distribution  systems.  With  the  majority  of  these  systems  installed 
more  than  25  years  ago,  an  extensive  repair  or  replacement  can  be  anticipated 
due  to  the  deteriorated  pipe  insulation  and  the  severely  corroded  carrier 
pipes  and  conduit  casings.  Nonmetallic  piping  materials  can  overcome  the 
corrosion  problems  experienced  with  steel  piping.  However,  most  nonmetallic 
materials  have  fairly  low  operating  temperature  limits  and  become  softer 
or  decomposition  at  elevated  temperatures.  In  direct  burial  conduit  heat 
distribution  systems,  hot  water  or  steam  are  conveyed  through  a pair  of 
parallel  buried  pipes,  which  are  covered  by  thermal  insulation  and  encased 
in  separate  exterior  conduits.  The  transmission  heat  losses  from  insulated 
piping  to  the  surrounding  earth  account  for  the  major  portion  of  the 
operating  costs. 
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Various  experimental  techniques,  such  as  field  testing  for  measuring  the 
heat  loss  from  the  underground  pipes,  are  expensive  and  time  consuming. 
Mathematical  modeling  is  an  alternative  approach,  which  provides  a relatively 
inexpensive  and  rapid  means  for  predicting  the  performance  of  heat 
distribution  system.  It  can  also  be  used  for  assessing  the  effects  of 
various  system  variables  such  as  pipe  size,  geometrical  configuration, 
insulation  thickness  and  operating  fluid  temperatures  on  the  system 
performance.  In  the  design  of  a new  underground  system  or  the  improvement 
of  the  existing  one,  mathematical  modeling  is  a valuable  tool  as  any 
modifications  to  the  design  can  be  implemented  in  the  model  and  tested 
easily.  The  predicted  results  from  mathematical  modeling,  for  example, 
the  surface  temperature  of  outer  conduit  casing,  can  be  used  as  a guide 
for  selecting  proper  conduit  and  coating  materials  for  the  underground 
systems.  The  prediction  of  the  heat  loss  from  the  buried  pipes  to  the 
surrounding  earth  is  information  needed  in  pipeline  design.  This  is  due 
to  heat  transfer  to  or  from  the  pipes  which  has  a significant  affect  to 
both  the  pumping  requirements  of  the  fluid  carried  through  the  pipes  and 
the  selection  of  optimum  pipe  insulation  thicknesses  based  on  life-cycle 
cost  analysis.  The  determination  of  heat  losses  and  surface  temperatures 
of  heat  distribution  systems  can  provide  a basis  for  refinement  of  guide 
specification  for  improved  guidance  on  design,  construction,  and  installation 
procedures  of  efficient  and  reliable  underground  systems. 

This  report  presents  the  computational  procedures  to  predict  the  heat 
losses  and  temperature  distributions  within  and  around  two  insulated  pipes 
encased  in  the  same  conduit  casing  and  in  separate  casings.  The  report 
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also  describes  the  theoretical  basis  of  computer  models  which  are  developed 
based  on  the  finite  element  technique  to  solve  a two-dimensional  steady- 
state  heat  transfer  between  the  surfaces  of  two  parallel  buried  pipes  and 
the  ground  surface.  A comparison  of  the  predicted  results  from  the 
computer  models  with  the  field  data  obtained  from  the  directly  buried 
conduit  heat  distribution  system  installed  at  Fort  Jackson,  SoutK  Carolina, 
by  the  U.  S.  Army  Cold  Regions  Research  and  Engineering  Laboratory  (CRREL) 
is  given  in  this  report. 

2.  Theoretical  Basis 

The  heat  transfer  to  the  surrounding  soil  from  two  insulated  pipes  installed 
horizontally  and  encased  in  a conduit  casing,  which  is  buried  at  a finite 
depth  into  the  earth,  involves  complex  geometrical  configurations,  composite 
materials  and  three  modes  of  energy  transport  including  conduction, 
convection  and  radiation.  With  this  nonlinear  heat  transfer  problem,  it 
is  not  possible  to  obtain  closed-form  analytical  solutions  for  the  rates 
of  heat  flow  in  the  vicinity  of  a directly  buried  conduit  system.  To  obtain 
approximate  solutions  to  the  differential  equations  governing  heat  flows, 
the  finite  element  method  is  used  because  of  its  adaptability  to  irregularly 
shaped  geometries  and  flexibility  for  different  element  sizes. 

2 . 1 Governing  Equation  and  Boundary  Conditions 

With  reference  to  the  Cartesian  coordinate  system,  the  differential 
equation  describing  two-dimensional  heat  conduction  within  an  isotropic 
solid  under  steady- state  conditions  with  no  internal  heat  generation  is 
given  by: 
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(1) 


with  the  constant  temperature  specified  on  a portion  of  the  boundary; 


T = Tb  ; on  $1 


(2) 


and  the  convection  boundary  conditions  specified  on  a part  of  the  boundary: 


(3) 


where  T is  temperature,  kj  is  thermal  conductivity  in  the  j -direction,  x 
and  y are  Cartesian  coordinates,  Ti,  is  the  prescribed  temperature  for  the 
boundary  segment  Sj^,  $2  is  the  boundary  segment  subject  to  convective  heat 
transfer,  Ix  and  ly  are  direction  cosines  of  a vector  perpendicular  to  S2. 
h is  the  convective  heat  transfer  coefficient,  and  is  the  temperature 
of  the  external  environment. 

2.2  Method  of  Solution 

Three -node  triangular  elements  are  used  along  with  a linear  variation  of 
temperatures  within  each  element  as  the  shape  function  to  discretize  the 
domain  of  interest.  Since  the  temperature  distribution  Is  a function  of 
space,  the  local  temperature  within  each  triangular  element  can  be 
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expressed  of  the  form: 


T(x,y)  = [N(x.y)]  (Ti)  (4) 

where  N(x,y)  and  are  the  shape  fxmctions  and  temperature  at  node  i, 
respectively. 

Using  the  approximation  for  the  unknown  temperatures  as  given  by  equation 
4 and  the  Galerkin  weighted  residual  method  [1-4],  the  heat  conduction 
equation  can  be  changed  into  a system  of  simultaneoixs  equations,  which  can 
be  written  in  a matrix  form  as: 

(K)  (T)  = (F)  (5) 

where  [K]  is  the  conductance  matrix,  (T)  is  the  column  vector  of  unknown 
nodal  temperatures,  and  (F)  is  the  column  matrix  consisting  of  the  boundary 
conditions . 


The  typical  elements  of  the  matrix  and  vectors  in  equation  5 can  be 
expressed  as : 


aN^  aN 
Kij  = Sv  ar  ^ 


aN.  aN 

+ k ^ ^ ) dV  + J 


y ay  ay 


S2  h Ni  Nj  d S 


(6) 


Fi  = /s2  h Ta  Ni  ds 


(7) 


Assemblage  of  these  element  equations  followed  by  modification  of  the 
assembled  equations  to  account  for  constant  temperature  and  convective  heat 
transfer  boundary  conditions  gives  the  global  system  of  equations.  This 
system  of  linear  equations  is  solved  for  the  unknown  temperature  vector  by 
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the  LU  decomposition  method  [5].  In  this  method,  the  conductance  matrix 
Is  factored  Into  the  product  of  a lower- triangular  matrix  and  an  upper- 
triangular  matrix.  The  resulting  set  of  equations  Is  then  solved  using 
forward  substitution  followed  by  backward  substitution. 

The  monthly  average  earth  temperature  for  a given  depth  used  as  the 
prescribed  temperature  along  the  outermost  perimeter  of  the  earth  region 
Is  a function  of  the  site  location,  month  of  year,  and  the  depth  below  the 
ground  surface.  It  can  be  estimated  from  the  following  equation  consisting 
of  a harmonic  function  [6] : 


where  T = the  monthly  average  earth  temperature,  °C  (°F) 

T^  = the  annual  average  earth  temperature  of  the  site,  °C  (°F) 

T^j  = the  annual  amplitude  of  the  monthly  average  temperature  cycle, 
OC  (F) 

y = the  depth  from  the  ground  surface,  m (ft) 
w = angular  frequency  of  the  annual  cycle,  rad/h. 
a = thermal  diffusitivity  of  the  soil,  m^/j^  (ft^/h) 
t = the  elapsed  time  from  January,  month 

Heat  transfer  by  free  convection  in  the  airspace  and  radiative  exchange 
between  the  surfaces  bounded  by  the  Insulated  pipe  and  its  horizontal 
conduit  casing  is  treated  as  a fictitious  non-transparent  material  having 
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equivalent  combined  convective  and  radiative  conductances  in  energy 
transport  process.  The  rate  of  convective  transfer  can  be  expressed  in 
terms  of  an  equivalent  heat  conduction  of  the  form  [7,8]; 

q = Ke  (Th  - Tc)/L  (9) 

where  q = the  average  heat  flow'  rate,  W/m^  (Btu/h*ft^^ 

Kg  = the  effective  thermal  conductivity  of  the  enclosed  air 
layer,  W/m.K  (Btu/h* ft . °F) 

= the  temperature  of  the  hot  surface,  °C(°F) 

Tg  = the  temperature  of  the  cold  surface,  qC(qF) 

L = the  characteristic  thickness  of  air  layer,  m(ft) . 

An  effective  thermal  conductivity,  which  is  defined  as  the  thermal 
conductivity  of  stationary  air  in  an  enclosed  space  necessary  to  transfer 
the  same  amount  of  heat  as  the  moving  air,  is  used  to  account  for  free 
convection  in  air  confined  between  the  heated  insulated  pipe  and  the 
cooler  conduit  casing.  The  ratio  of  the  effective  to  actual  thermal 
conductivity  of  the  enclosed  airspace  is  a complex  function  of  geometrical 
configuration  of  wall  surfaces,  the  thermal  properties  and  motion  of 
confined  air,  the  characteristic  dimension  or  the  thickness  of  the  air 
layer,  and  the  temperature  difference  between  the  hot  and  cold  surfaces  [8]. 

Free  convection  heat  transfer  in  the  annular  space  between  two  long 
horizontal  concentric  cylinders  was  studied  experimentally  using  optical 
methods  and  theoretically  using  numerical  techniques  by  Kuehn  and  Goldstein 
[9].  Their  experimental  results  using  air  and  water  as  the  working  media 
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were  correlated  in  terms  of  the  effective  thermal  conductivity.  The 
average  effective  to  actual  thermal  conductivity  for  air  can  be  expressed 
as  below: 

Ke/Kf  = 0.159  Ra  0.272  ^0) 

where  Kf  = thermal  conductivity  of  air  at  Tf 
Ra  = Rayleigh  number,  PrGr 
Pr  = Prandtl  number,  v/a 
u = Kinematic  viscosity  of  air 
a = thermal  diffusivity  of  air 
Gr  = Grashof  number,  g/SCTi-Tj,)  L^/u^ 
g = acceleration  of  gravity 
^ = thermal  expansion  coefficient  at  Tf 

L = gap  width,  L = Ro"Ri 

Ri,Ro  = radius  of  the  inner  and  outer  cylinder,  respectively 
Tj^,TQ  = temperature  of  the  inner  and  outer  cylinder,  respectively 

For  the  case  of  two  long,  horizontal  pipes  encased  in  the  same  conduit 
casing,  no  specific  experimental  data  or  correlation  results  on  free 
convection  in  confined  space  are  available.  In  the  absence  of  experimental 
correlation,  the  heat  transfer  by  free  convection  in  airspaces  between  two 
small  eccentric  horizontal  cylinders  and  a large  cylinder  may  be  treated 
similarly  as  natural  convection  around  a horizontal  cylinder.  It  can  be 
estimated  by  an  expression  of  the  form  [10]. 

0.25 

Nu  = 0.53  RaD  (H) 

where  Nu  = Nusselt  number,  hL/Kf=  Ke/Kf 
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h = average  heat  transfer  coefficient 
L = hypothetical  gap  width,  L = (Dq  - D)/2 
D = effective  diameter  of  the  two  small  cylinders, 
which  is  defined  as  the  diameter  of  a cylinder 
having  the  perimeter  equal  to  the  total  perimeter 
of  the  small  cylinders. 

Dq  = the  diameter  of  the  large  cylinder 
RaD  = Rayleigh  number,  Pr  GrD 
Crj)  = Grashof  number,  g^CT^-To) 

The  radiant  exchange  between  the  outer  surfaces  of  the  insulated  pipes  and 
the  inner  surface  of  the  conduit  casing  is  treated  as  a multiple  radiative 
interchange  between  two  infinitely  long  concentric  cylinders  with  diffuse  - 
gray  surfaces  of  different  emissivities  and  temperatures.  The  equivalent 
thermal  conductance  due  to  radiative  transfer  to  be  added  to  the  convective 
component  to  yield  the  overall  conductance  can  be  determined  from  the 
following  expression: 

Kj.  = hj-L  (12) 

where  hj,  = a (Ti  + Tq)  (Ti2  + Tq^)/  [l/e^  + (Di/Do)  (Dq)  (l/co  " D] 
a = the  Stefan-Boltzmann  constant 
£i,eo  = bhe  total  emissivity  of  the  exposed  surface  of 
the  inner  and  outer  cylinder,  respectively 
Di,DQ  = the  outside  diameter  of  the  inner  cylinder  and  the 
inside  diameter  of  the  outer  cylinder,  respectively 
L = the  gap  width,  (DQ-Dj^)/2 

For  an  insulated  piping  system,  the  surface  film  resistance  between  the  hot 
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fluid  and  the  steel  pipe  and  the  thermal  resistance  of  pipe  wall  can 
generally  be  neglected  compared  to  the  thermal  resistance  of  the  insulation. 
The  rate  of  heat  loss  from  an  insulated  pipe  can  be  calculated  from  the 
following  equation  along  with  the  calculated  value  of  average  temperature 
drop  across  the  pipe  insulation  layer: 

q = 2^rki(Ti  - T2)/ln  (r2/ri)  (13.) 

where  q = the  heat  loss  rate  per  unit  length  of  the  insulated  pipe, 

W/m  (Btu/h.ft) 


^i  “ 

^1.^2  = 

Tl 


T2 


thermal  conductivity  of  insulation  material,  W/m.K 


(Btu.in/h.ft.2  Op) 
the  inside  and  outside  radius 
respectively,  m(ft) 
the  surface  temperature  of  the 
radius  , which  is  assumed  to 
fluid  temperature, 
the  surface  temperature  of  the 
radius,  °C  (°F) 


if  the  insulation  layer, 

insulation  layer  at  inner 
be  the  same  as  the  working 

insulation  layer  at  outer 


The  thermal  conductivity  of  pipe  insulation  is  generally  a function  of  its 
mean  temperature.  For  calcium  silicate  insulation,  the  thermal  conductivity 
function  is  based  on  a look-up  table  and  stored  in  a subprogram  to  provide 
the  temperature  dependent  thermal  conductivity  value.  The  thermal 
conductivity  of  mineral  wool  insulation  is  obtained  using  a least-squares 
regression  analysis  of  the  test  data  from  the  pipe  insulation  used  for 
the  underground  systems  installed  at  Fort  Jackson,  South  Carolina.  It 
can  be  expressed  by 
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(14) 


Ki  = 0.242  + 1.5501  x 10  Tn,  + 7 . 500  x 10  Tm^ 
where  = thermal  conductivity  of  the  insulation,  Btu-in/h.ft^  °F 

Tju  = the  mean  temperature  of  the  insulation,  °F 


The  thermal  conductivity,  kinematic  viscosity  and  Prandtl  number  of  air 
are  a function  of  its  mean  temperature  and  can  respectively  be  approximated 
as 


Ka  = 0.01319  + 2.50  x 10 

ua  = 1.2624  X 10  + 5.40  X 10  Ta 

Pr  = 0.7185  - 1.275  x 10  Ta 


(15) 


where 

Ka  = thermal  conductivity  of  air,  Btu/h.ft.°^ 

Ta  = the  bulk  air  temperature  within  the  airspace,  °F 
Ua  = kinematic  viscosity  of  air,  ft^/s 


The  bulk  air  temperature  can  be  approximated  by  the  average  of  the  outer 
surface  temperature  of  the  insulation  and  the  inner  surface  temperature  of 
the  conduit  casing.  For  the  case  of  two  pipes  encased  in  the  same  conduit, 

an  effective  insulation  outer  surface  temperature  is  defined  to  be  equal 
to  the  weighted  average  of  the  external  surface  temperatures  of  two  insulated 
pipes . 


3.  Description  of  Computer  Programs 

Two  finite  element  computer  programs  called  DIRECTl  and  DIRECT2  have  been 
developed  for  numerical  modeling  of  two-dimensional,  steady-state  heat 
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conduction.  This  involves  directly  buried  conduit  distribution  systems 
with  two  parallel  pipes  encased  in  the  same  conduit  casing  and  in  separate 
casings,  respectively.  Each  computer  program  consists  of  a main  program 
and  nine  subroutines  and  the  computer  codes  are  written  in  FORTRAN 
language.  In  addition  to  reading  in  portions  of  input  data,  the  main 
program  calls  pertinent  subroutines,  calculates  the  conductance  matrix 
and  excitation  vector  modified  for  the  given  convection  and  constant 
temperature  boundary  conditions,  and  prints  out  the  calculated  nodal 
temperatures . 

Subroutine  TGO  calculates  the  average  undisturbed  earth  temperatures  at 
various  depths  for  the  month  of  interest.  Subroutine  PIPEO  is  used  to  read 
in  the  dimensions  of  the  insulated  pipes  and  the  inner  and  outer  regions 
of  the  surrounding  earth,  echos  the  data  to  allow  the  checking  of  input  data, 
and  generates  the  rectangular  coordinates  for  each  nodal  point  of  the  two- 
pipe  system.  Subroutines  INSULK  and  SOILK  provide  the  thermal  conductivity 
of  pipe  insulation  including  calcium  silicate  and  mineral  wool,  and  soil 
thermal  conductivity,  respectively  by  linear  interpolation  of  sets  of 
thermal  conductivity  versus  mean  temperature  data.  Subroutines  TGXX 
furnish  the  boundary  temperatures  of  the  outer  earth  region  and  TWOPIP 
calculates  the  rate  of  heat  loss  from  the  two  insulated  pipes  to  the  soil 
surrounding  the  heat  distribution  system.  Subprograms  EQUIKO  determines 
the  equivalent  thermal  conductivity  of  airspace  within  the  pipe  and  SOLVE 
is  used  to  solve  system  of  simultaneous  equations  by  LU  decomposition 
method  [5].  Subroutine  PIPEHL  calculates  the  average  temperature  drops 
across  the  pipe  insulation  layers,  the  rates  of  heat  loss  from  both 
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underground  pipes,  and  prints  out  the  results  of  these  calculations. 

A computer  program  called  DBJACKS  similar  to  the  computer  code  DIRECTl  in 
program  structures  and  flow  control,  has  been  developed  to  simulate  the 
conditions  of  thermal  measurements  performed  by  the  U.S.  Army  Cold  Regions 
Research  and  Engineering  Laboratory  for  the  high  temperature  hot  water 
distribution  system  installed  at  Fort  Jackson,  South  Carolina.  A listing 
of  the  source  codes  of  these  computer  programs  is  given  in  Appendix  B. 

In  using  the  computer  codes,  the  underground  systems  to  be  modeled  numerically 
are  first  divided  into  regions  of  interest.  These  include  pipe  insulation 
covering  the  carrier  pipes,  airspace  within  the  pipes,  the  conduit  casing, 
back-fill  region,  and  the  inner  and  outer  earth  regions  in  the  vicinity  of 
the  distribution  system.  These  regions  are  then  discretized  into  triangular 
elements.  To  increase  the  degree  of  calculation  accuracy,  the  smaller 
element  sizes  are  used  in  the  areas  of  anticipated  higher  temperature 
gradients.  The  larger  elements  are  employed  for  the  areas  having  smaller 
temperature  gradients.  In  labeling  of  the  nodal  points,  the  nodes  along 
the  inner  and  outer  boundaries  with  the  specified  temperatures  and  convective 
heat  fluxes  are  numbered.  This  is  done  first  to  obtain  a reduced  number 
of  simultaneous  equations,  which  will  be  solved  for  the  unknown  nodal 
temperatures . 

The  finite  element  meshes  for  the  direct  burial  conduit  systems  for 
distributions  of  steam  and  hot  water  are  shown  in  Figures  1 through  3. 

These  figures  are  divided  into  three  sets:  the  "a"  figures  illustrate  the 
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underground  system  consisting  of  both  pipes  encased  in  separate  conduits 
the  "b"  figures  show  the  distribution  system  involving  two  pipes  covered 
by  the  same  outer  conduit,  and  the  "c"  figures  illustrate  the  underground 
system  installed  at  Fort  Jackson  and  having  both  pipes  encased  in  the  same 
conduit.  Figures  l.a  to  l.c  show  the  finite  element  design  for  the  inner 
and  outer  earth  regions  surrounding  the  distribution  pipes.  The  finite 
element  grids  for  back-fill  region  adjacent  to  the  underground  system  are 
shown  in  Figures  2. a through  2.c.  Figures  3. a through  3.c  illustrate  the 
finite  element  meshes  for  pipe  insulation,  the  conduit  casing,  and  airspace 
between  the  insulated  pipes  and  the  casing.  The  finite  element  meshes 
contained  in  the  computer  programs  consist  of  248  triangular  elements  and 
142  nodal  points  for  program  DIRECT2,  212  elements  and  124  nodes  for 
DIRECTl , and  224  elements  and  131  nodes  for  DBJACKS , respectively.  Table 
1 presents  a summary  of  finite  element  meshes  representing  regions  of  the 
major  components  of  underground  systems  where  both  pipes  are  encased  in 
separate  conduits  and  in  the  same  conduit. 

4.  Comparison  of  Calculated  Results  with  Experimental  Data 
The  direct  burial  conduit  heat  distribution  system  installed  at  Fort 
Jackson,  South  Carolina  and  modeled  numerically,  consists  of  two  5-inch 
(127  mm)  high  temperature  hot  water  pipes  with  the  heat  supply  pipe  placed 
directly  above  the  heat  return  pipe.  Each  pipe  is  insulated  by  a 1.5 
inches  (38  mm)  thick  layer  of  mineral  wool.  The  hot  water  supply  and 
return  lines  are  encased  in  a 20  inches  (508  mm)  outside  diameter  steel 
conduit  having  a wall  thickness  of  1/8  inch  (3.2  mm).  Both  pipes  are 
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separated  by  a distance  of  9.5  inches  (235  mm)  between  pipe  centers.  The 
conduit  casing  is  located  at  3 ft  8 inches  (1.12  m)  beneath  the  ground 
surface.  In  numerical  simulation,  the  piping  system  is  assumed  to  be 
surrounded  by  earth  having  a thermal  conductivity  of  7.0  Btu.in/h.ft^  °F 
and  an  average  temperature  of  50  °F.  The  ground  surface  is  exposed  to 
ambient  air  of  seasonally  varying  temperatures.  All  experimental  data 
obtained  by  the  CRREL  from  the  test  site  were  recorded  at  6 - hour 
intervals.  Numerical  calculations  were  made  based  on  the  average  values 
over  a two-day  period  of  the  measured  heat  supply  and  return  pipe 
temperatures,  and  the  experimentally  determined  undisturbed  earth 
temperatures  at  various  depths . 

Two  input  data  files  SDTAIFJ  and  SDTA2FJ  are  created  prior  to  execution  of 
computer  program  DBJACKS  for  numerical  simulations  of  the  thermal  performance 
of  the  directly  buried  conduit  heat  distribution  system  installed  at  Fort 
Jackson.  These  input  data  files  along  with  the  outputs  from  the  computer 
program  are  presented  in  Appendix  A.  The  SDTALFJ  shown  in  Appendix  A.l 
refers  to  the  main  program  and  subroutines  PIPEJ  and  TGO  for  data  input. 

It  contains  the  title  of  the  computer  rtin,  total  numbers  of  nodal  points 
and  triangular  elements.  The  first  element  number  of  pipe  insulation,  the 
first  node  numbers  of  pipe  insulation  for  the  heat  supply  and  return  pipes, 
and  the  month  of  interest  are  also  included.  In  addition  data  for  run 
control  parameters,  thermal  conductivity  and  dimensions  of  the  inner  earth 
region,  estimated  average  air  temperature,  and  temperature  drop  across 
airspace  are  presented.  Also  shown  are  thermal  conductivities  of  the 
conduit  casing  and  the  back-fill  soil,  the  total  emissivities  and  average 
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temperatures  of  outer  surface  of  pipe  insulation  and  inner  surface  of  the 
conduit  casing,  the  pipe  fluid  temperatures,  the  material  type,  thermal 
conductivity  and  thicknesses  of  pipe  insulation,  the  pipe  sizes  and 
locations  of  the  pipe  centers,  the  size  and  depth  of  the  outer  casing, 
their  thermal  properties  and  dimensions  of  the  outer  earth  region,  and  the 
annual  average  earth  temperature  and  amplitude  of  the  monthly  temperature 
cycle  for  the  site  involved.  As  shovm  in  Appendix  A. 2,  the  element  data 
file  SDTA2FJ  consists  of  the  element  number,  the  node  numbers  for  its  three 
vertices  and  the  material  type  for  each  triangular  element,  total  number 
of  elements  subject  to  convection  boundary,  the  surface  convection 
coefficients  and  ambient  temperatures  for  three  sides  of  each  element 
that  experiences  convection  loss. 

The  boundary  conditions  for  numerical  simulations  are  the  constant 
temperatures  of  working  fluids  flowing  inside  the  carrier  pipes,  the 
specified  undisturbed  earth  temperatures  along  the  outermost  perimeter  of 
the  earth  region,  and  rate  of  convective  heat  transfer  between  the  ground 
surface  in  the  vicinity  of  the  underground  system  and  the  ambient  air. 

The  undisturbed  earth  temperatures  at  various  depths  for  a specific  month 
were  calculated  from  equation  8 in  which  the  two  parameters:  the  average 
earth  temperature  of  the  site  and  the  best- fit  of  the  measured  earth 
temperature  data  to  the  equation  through  the  method  of  least  squares. 

The  system  of  matrix  equations  is  solved  for  the  unknown  nodal  temperatures 
using  the  LU  decomposition  technique.  With  adjustment  of  insulation  and 
soil  thermal  conductivities  to  account  for  temperature  effect,  an 
iterative  procedure  is  used  until  the  heat  loss  from  the  underground 
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pipes  reaches  a steady-state  condition. 


Computer  calculations  of  pipe  heat  loss  and  temperature  distributions  in 
the  vicinity  of  the  underground  system  were  performed  on  the  mainframe 
Cyber  855  at  NIST.  These  calculations  required  approximately  5 to  8 
seconds  of  execution  time.  The  outputs  from  the  computer  program  included: 
the  rates  of  heat  loss  from  carrier  pipes,  the  resultant  temperature  at 
each  nodal  point,  average  temperature  drops  across  pipe  insulation  layers, 
equivalent  thermal  conductivities  of  airspace  due  to  free  convection  and 
radiative  transfer,  and  cartesian  coordinates  of  all  nodal  points.  A 
listing  of  output  file  SOUTFJ  from  program  DIRECFJ  for  the  sample  case  is 
given  in  Appendix  A. 3. 

Numerical  predictions  of  the  thermal  performance  of  this  high  temperature 
hot  water  distribution  system  were  carried  out  to  simulate  the  tests 
conducted  during  the  January  1-10,  1987  period.  A comparison  of  the 
predicted  temperature  distributions  within  and  around  the  two  insulated 
pipes  encased  in  the  same  metallic  conduit  with  the  corresponding  measured 
values  obtained  by  the  CRREL  at  the  test  site  is  tabulated  in  Table  2.  The 
average  value  and  the  range  of  the  calculated  surface  temperatures  at  8 
nodal  points  distributed  evenly  around  the  outer  surface  of  pipe  insulation, 
or  the  inner  surface  of  the  conduit  casing  are  also  presented  in  Table  2 
for  comparison.  In  general,  the  predicted  results  are  in  reasonably  good 
agreement  with  the  experimental  data.  The  difference  between  the  calculated 
and  measured  temperatures  for  all  measuring  locations  with  the  exception 
of  the  outermost  boundary  is  found  to  range  from  0 to  23  percent  with  an 
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average  of  7 percent.  A large  discrepancy  between  the  predicted  and 
experimental  temperatures  on  the  outer  surface  of  pipe  insulation  for  the 
heat  supply  line  can  be  explained  by  the  mean  equivalent  thermal  conductivity 
used  in  calculations.  It  is  too  large  compared  to  the  actual  local 
equivalent  conductivity  which  would  result  in  relatively  low  insulation 
surface  temperature.  Kuehn  and  Goldstein  [9]  reported  that  the  local 
equivalent  conductivity  or  local  heat  transfer  coefficient  of  horizontal 
concentric  cylinders  varied  with  angular  position.  The  highest  value 
occurred  at  the  bottom  of  the  inner  cylinder  with  the  smallest  value  at  the 
top.  With  the  absence  of  an  expression  for  the  local  heat  transfer 
coefficient,  it  is  relevant  to  use  the  average  value  of  the  predicted 
temperatures  on  a cylindrical  surface  in  correspondence  to  the  mean  heat 
transfer  coefficient  employed  in  numerical  calculations. 

As  illustrated  in  Table  2,  the  deviations  of  the  calculated  results  based 
on  average  surface  temperatures  from  the  measured  values  are  found  to  vary 
from  0 to  17  percent  with  a mean  of  6 percent.  Figure  4 shows  the  temperature 
data  predicted  with  the  finite  element  models  for  temperatures  at  selected 
locations  within  and  around  the  directly  buried  conduit  system  at  Fort 
Jackson  4 against  the  temperature  values  determined  by  thermocouples  at  the 
test  site.  A similar  plot  for  comparing  the  temperature  results  from 
thermocouple  measurements  with  the  predicted  values  including  the  average 
values  of  eight  local  surface  temperatures  on  the  outer  surface  of  the 
insulated  pipe  and  on  the  inner  surface  of  the  conduit  casing  is  shown  in 
Figure  5.  These  plots  show  that  a reasonably  good  correlation  exists 
between  the  temperature  values  obtained  from  the  model  predictions  and  the 
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field  measurements  because  all  the  data  points  lie  close  to  the  line  of 
perfect  agreement. 

5.  Comparisons  of  Finite  Element  and  Thermal  Analysis  Programs  with 
Guide  Specification  Table 

A thermal  analysis  computer  program  called  DIRECT?,  was  developed  based  on 
steady-state,  one -dimensional  heat  conduction  for  an  insulated  pipe  with 
an  airspace  conduit  buried  horizontally  beneath  the  ground  surface.  The 
theoretical  treatments  considered  the  condition  that  hot  fluid  flows  inside 
a steel  pipe  and  heat  is  transferred  radially  by  conduction  and  convection 
through  a series  of  cylindrical  shells.  These  consist  of  steel  pipe, 
thermal  insulation,  airspace  and  conduit  casing,  and  the  soil  layer  around 
the  underground  pipe.  The  calculation  procedure  is  essentially  similar  to 
that  used  to  derive  the  pipe  insulation  tables  for  the  maximum  permissible 
heat  loss  values  listed  in  Guide  Specification  CEGS-15705  [11].  A listing 
of  the  source  code  of  thermal  analysis  computer  program  is  given  in 
Appendix  B.  The  outputs  from  thermal  analysis  program  include  the  heat 
loss  from  a single  pipe  of  directly  buried  conduit  system,  the  mean 
insulation  thermal  conductivity,  and  the  inner  and  outer  surface  temperatures 
of  both  the  pipe  insulation  and  the  conduit  casing. 

Some  calculations  using  the  finite  element  computer  simulation  program, 
DIRECT2,  and  the  thermal  analysis  program  were  carried  out  on  t3rpical 
direct  burial  systems  with  nominal  6-in.  (152  mm)  and  3-in.  (76  mm)  pipes 
and  different  working  fluid  temperatures  ranging  from  200  °F  to  400  °F. 

The  results  on  pipe  heat  loss  were  compared  with  the  corresponding  maximum 
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permissible  heat  loss  values  given  in  Table  5 of  Guide  Specification  CEGS- 
15705.  In  these  calculations,  the  heat  supply  and  return  pipes  were 
located  side  by  side  at  4.0  ft.  (1.22m)  below  the  ground  surface,  and 
separated  by  a distance  of  1.5  ft.  (0.46m)  between  pipe  centers.  The 
thicknesses  of  calcium  silicate  insulation  used  for  the  pipes  were  varied 
from  1.5  inches  (38mm)  to  3.0  inches  (76  mm).  An  average  ground  temperature 
of  55  °F  (12. 8C)  and  the  thermal  conductivity  of  soil  equal  to  5 
Btu. in/h. ft^ . °F  (0.72  W/m.C)  were  used  in  calculations.  The  thicknesses 
of  outer  casing  and  airspace  between  the  insulated  pipe  and  the  casing 
were  1/8  inch  (3mm)  and  about  1 inch  (25mm),  respectively. 

The  pipe  heat  losses  predicted  by  the  finite  element  and  thermal  analysis 
computer  programs  for  directly  buried  systems  for  hot  water  and  steam 
distributions  are  shown  in  Table  3.  In  general,  for  the  heat  supply  pipes, 
the  calculated  heat  loss  values  from  the  finite  element  model  are 
approximately  5 percent  smaller  compared  to  the  values  listed  in  the 
table  of  the  guide  specification.  For  the  heat  return  pipes  having  the 
identical  pipe  size  as  the  heat  supply  pipes,  the  finite  element  model 
gives  about  22  percent  smaller  heat  loss  values  in  comparison  to  the 
guide  specification  table.  For  the  case  of  steam  distribution  as 
illustrated  in  the  table,  this  difference  in  heat  loss  values  for  the 
heat  return  pipe  become  greater  with  the  increased  size  of  the  heat 
supply  pipe.  The  rate  of  heat  loss  from  a single  pipe,  especially  for 
the  heat  return  pipe,  is  considerably  influenced  by  the  pipe  size  and 
working  fluid  temperature  of  its  neighboring  pipe,  which  serves  as  a heat 
source  or  heat  sink  to  continuously  add  or  remove  a portion  of  heat  from 
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the  thermal  field  involved.  The  heat  loss  results  calculated  by  thermal 

analysis  program  are  consistently  about  3 percent  greater  than  the 

maximum  allowable  heat  loss  values  tabulated  in  the  guide  specification  table. 

6.  Conclusions 

The  temperature  distributions  in  underground  heat  distribution  systems, 
which  consist  of  a pair  of  insulated  pipes  encased  in  the  same  conduit  and 
in  separate  conduits  within  the  surrounding  soil,  and  the  heat  loss  rates 
from  the  buried  pipes  were  calculated  using  computer  simulation  programs 
utilizing  the  finite-element  method. 

General  formulation  of  the  relevant  equations  governing  heat  flow  and 
boundary  conditions  for  steady-state,  two-dimensional  heat  conduction 
problems  are  presented.  The  computational  scheme,  the  input  data  required 
for  executing  the  simulation  programs,  and  the  outputs  from  the  computer 
iruns  for  sample  cases  are  described.  The  spatial  distributions  of 
temperatures  predicted  by  the  finite  element  models  for  a high  temperature 
hot  water  distribution  system  installed  at  Fort  Jackson  involving  two 
insulated  pipes  encased  in  the  same  conduit  are  generally  consistent  with 
the  experimental  data  obtained  by  the  CRREL  at  the  test  site.  The 
difference  between  the  predicted  and  measured  temperatures  is  found  to 
range  from  0 to  17  percent  with  an  average  of  approximately  6 percent. 

A comparison  was  made  of  the  pipe  heat  losses  predicted  by  the  finite 
element  and  the  thermal  analysis  models  for  a range  of  pipe  working  fluid 
temperatures  and  insulation  thicknesses  with  the  maximum  permissible  heat 
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loss  values  listed  in  Table  5 of  Guide  Specification  CEGS- 15705.  The 
finite  element  model  gives  approximately  5 percent  smaller  heat  loss  value 
for  the  heat  supply  pipe  and  22  to  27  percent  less  in  heat  loss  for  the 
heat  return  pipe  compared  to  the  guide  specification  table,  depending  upon 
the  pipe  size,  separation  distance  and  working  fluid  temperature  of  the 
neighboring  pipe.  It  is  anticipated  that  the  pipe  located  in  the  vicinity 
serves  as  a heat  source  or  heat  sink  to  interact  with  the  thermal  field 
involved.  The  pipe  heat  loss  values  predicted  by  the  thermal  analysis 
model  derived  based  on  steady-state,  one -dimensional  radial  heat  conduction 
through  a single  insulated  pipe  are  found  to  agree  within  3 percent  of  the 
values  tabulated  in  the  guide  specification.  Experimental  data  and 
correlation  equations  of  local  heat  transfer  coefficient  are  needed  for  a 
more  accurate  prediction  of  heat  flows  in  an  airspace  between  the  insulated 
pipe  and  the  inner  surface  of  the  conduit  casing. 
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Table  1 A Summary  of  Finite  Element  Meshes  for  Underground  Systems  with  Two  Pipes  Encased 
in  Separate  Conduits  and  in  the  Same  Conduit 


Pipes  in  Separate 
Conduits 

Pipes  in  Same 
Conduits 

Simulation 

Hot  Water 

System 

1 . Inner  Earth 

Region 

Elements  1 through 

32  (Fig.  l.a) 

Elements  1 through 

32  (Fig.  l.b) 

Elements  1 through 

32  (Fig.  l.c) 

2 . Outer  Earth 

Region 

Elements  33  through 

90  (Fig.  l.a) 

Elements  33  through 

90  (Fig.  l.b) 

Elements  33  through 
96  (Fig.  l.c) 

3.  Pipe 

Insulation 

Elements  91  through 

122  (Fig.  3. a) 

Elements  91  through 
122  (Fig.  3.C) 

Elements  97  through 
128  (Fig.  3.C) 

4.  Airspace 

Between  Insulated 
Pipe  and  Conduit 
Casing 

Elements  123  through 
154  (Fig.  3. a) 

Elements  123  through 
148  (Fig.  3.b) 

Elements  129  through 
160  (Fig.  3.C) 

5.  Conduit 

Casing 

Elements  155  through 
186  (Fig.  3. a) 

Elements  149  through 
164  (Fig.  3.b) 

Elements  161  through 
176  (Fig.  3.C) 

6.  Backfill 

Region 

Elements  187  through 
248  (Fig.  2. a) 

Elements  165  through 
212  (Fig.  2.b) 

Elements  177  through 
224  (Fig.  2.C) 

Total  Number  of 

Nodal  Points: 

142 

124 

131 

Total  Number  of 
Elements : 

248 

212 

224 
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TABLE  2.  COMPARISON  OF  CALCULATED  RESULTS  WITH  EXPERIMENTAL  DATA 


OUTDOOR  SURFACE  TEMP.  (DEG.  AIRSPACE 


AIR 

TOP 

BOTTOM 

CONDUIT 

TEMP. 

(F) 

DATE 

VALUE 

TEMP.  rF') 

PIPE 

INSUL. 

PIPE 

INSUL. 

(INSIDE) 

TOP 

MID 

BOT 

1/1- 

Meas . 

41 

325 

156 

222 

127 

121 

140 

146 

116 

1/2/ 

Calc . 

122(131) 

134(127) 

118(117), 

120 

148 

117 

87 

[122-148] 

[118-147] 

[115-118] 

Diff.  (%) 

21.8(16.0) 

5. 5(0.0) 

2. 5(3. 3) 

14.3 

1.4 

0.9 

1/3- 

Meas . 

39 

323 

153 

206 

122 

119 

140 

142 

113 

1/4/ 

Calc . 

118(127) 

129(123) 

115(113) 

117 

143 

113 

87 

[118-143] 

[114-142] 

[111-115] 

Diff.  (%) 

22.9(17.0) 

5. 7(0. 8) 

3. 4(5.0) 

16.4 

0.7 

0.0 

1/5- 

Meas . 

42 

322 

152 

213 

122 

118 

138 

142 

113 

1/6/ 

Calc . 

120(129) 

131(125) 

115(113) 

119 

145 

115 

87 

[120-145] 

[116-144] 

[111-115] 

Diff.  (%) 

21.1(15.1) 

7. 4(2. 5) 

2. 5(4. 2) 

13.8 

2.1 

1.8 

1/7- 

Meas . 

46 

323 

154 

220 

124 

118 

139 

144 

114 

1/8/ 

Calc . 

123(131) 

134(127) 

119(117) 

121 

148 

118 

87 

[123-148] 

[118-147] 

[115-119] 

Diff.  (%) 

20.1(14.9) 

8. 1(2.4) 

0.8(0. 8) 

12.9 

2.8 

3.5 

1/9- 

Meas . 

48 

321 

153 

225 

125 

118 

139 

144 

114 

1/10 

Calc. 

124(133) 

135 

120(119) 

122 

149 

119 

/87 

[124-150] 

[120-149] 

[117-120] 

Diff.  (%) 

19.0(13.1) 

8. 0(3. 2) 

1.7(0. 8) 

12.2 

3.5 

4.4 

Note:  The  numerical  values  inside  the  parentheses  and  the  brackets 

are  the  average  value,  and  the  range  of  the  calculated 
surface  temperatures  at  8 nodai  points  distributed  evenly 
around  the  cylindrical  surface,  respectively. 
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EARTH  UNDISTURBED 


TEMP 


DATE 

VALUE 

2-in 

DEEP 

1/1- 

Meas . 

56 

1/2/ 

Calc . 

57 

87 

Diff. 

(%) 

1.8 

1/3- 

Meas . 

56 

1/4/ 

Calc . 

54 

87 

Diff. 

(%) 

3.6 

1/5- 

Meas . 

56 

1/6/ 

Calc . 

57 

87 

Diff. 

(%) 

1.8 

1/7- 

Meas . 

56 

1/8/ 

Calc . 

60 

87 

Diff. 

(%) 

7.1 

1/9- 

Meas . 

56 

1/10 

Calc . 

62 

/87 

Diff. 

(%) 

10.7 

(F) 

EARTH 

TEMP 

9 -in. 

3 -in. 

26 -i: 

DEEP 

DEEP 

DEE 

66 

46 

50 

61 

47 

49 

7.5 

2.1 

2.0 

66 

45 

49 

59 

46 

48 

10.6 

2.2 

2.0 

66 

44 

48 

61 

45 

47 

7.6 

2.3 

2.1 

66 

45 

49 

65 

46 

47 

1.5 

2.2 

4.1 

66 

47 

50 

66 

48 

49 

0.0 

2.1 

2.0 

SYSTEM 

56 -in.  HEAT  LOSS 

DEEP  fBtu/h-ft') 

51 

51  110 

0.0 

50 

50  106 

0.0 

49 

49  107 

0.0 

49 

49  108 

0.0 

50 

51  108 

2.0 
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Table  3.  Comparison  of  Pipe  Heat  Loss  Values  (Btu/h.ft)  Predicted 
by  Finite  Element  Computer  Program  with  Calculated 
Values  Obtained  from  Thermal  Analysis  Program  and 
Guide  Specification  Table 

Pipe  Working  Fluid  Temperature  (Degree  F) 


Calculated 
by 


200 


250 


350 


400 


1,  3-in.  Pipe: 


FECP 

29.8 

36.2 

61.7^ 

60.7b 

TACP 

(25.5)a 

(22.9)e 

33.8 

(33.1)b 

(28.8)f 

39.9 

62.5 

61.7 

CEGS 

32.7 

38.7 

60.6 

59.5 

Insulation 

Thickness ( in. ) 

: 1.5 

2.0 

2.0 

3.0 

2.  6-in.  Pipe 

FECP 

42.1 

57.4 

80.0c 

93.9^ 

(34.3)c 

(49.3)d 

81.8® 

97. 5f 

TACP 

48.8 

67.8 

81.5 

97.2 

CEGS 

47.8 

65.2 

80.3 

95.8 

Insulation 

Thickness (in. ) 

: 1.5 

1.5 

2.5 

2.5 

Notes : 

1 . FECP 

and  TACP 

represent 

finite  element 

and  thermal 

computer  programs,  respectively.  CEGS  denotes  Guide  Specification 
CEGS- 15705,  which  contains  the  heat  loss  values  listed  in  Table  5. 


2.  The  heat  loss  values  in  parentheses  are  calculated  based  on  the 
pipe  used  as  the  heat  return  pipe. 

3.  Superscripts  a through  f denote  pairs  of  a two-pipe  system  in 
calculations . 
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Figure  l.a  Finite  Element  Mesh  for  Inner  and  Outer  Earth  Regions  with 
Both  Pipes  in  Separate  Conduits 
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Figure  l.b  Finite  Element  Mesh  for  Inner  and  Outer  Earth  Regions  with 
Both  Pipes  in  the  Same  Conduit 
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Figure  1,C  Finite  Element  Mesh  for  Inner  and  Outer  Earth  Regions  for 
Underground  System  Installed  at  Fort  Jackson 
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Figure  2. a Finite  Element  Mesh  for  Back-fill  Region  Surrounding  Insulated 
Pipes  Encased  in  Separate  Conduits 
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Figure  2.b  Finite  Element  Mesh  for  Back-fill  Region  Surrounding  Insulated 
Pipes  Encased  In  the  Same  Conduit 
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Figure  2.C  Finite  Element  Mesh  for  Back-fill  Region  Around  Underground 
System  Installed  at  Fort  Jackson 
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Figure  3. a Finite  Element  Design  for  Pipe  Insulation,  Airspace  Layer  and 
Conduit  Casing  for  Two  Pipes  in  Separate  Conduits 
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Figure  3.b  Finite  Element  Design  for  Pipe  Insulation,  Airspace  Layer  and 
Conduit  Casing  for  Two  Pipes  in  the  Same  Conduit 
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Figure  3.C  Finite  Element  Design  for  Pipe  Insulation,  Airspace  Layer  and 
Conduit  Casing  for  Underground  System  Installed  at  Fort  Jackson 
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Figure  4.  A Comparison  Between  the  Temperature  Values  Predicted  with 
Finite  Element  Models  and  the  Values  Obtained  from  Thermocouple 


200 


Cd  &Jn^Dj»dujdj,  pejn sDd^ 


39 


Calculated  Temperature  (Deg  F) 

Figure  5.  A Comparison  of  the  Calculated  Temperatures  Including  the 
Average  Values  for  Cylindrical  Surfaces  with  the  Measured 


APPENDIX  A.  The  Input  Data  Files  and  the  Outputs  from  the  Computer  Programs 


A.l  A Listing  of  SDTAIFJ  Input  File  for  Program  DBJACKS 

HOT  WATER  SYSTEM(PIPES  IN  ONE  CASINC)AT  FORT  JACKSON.  1/7/87 

131 ,224,31 .97. 85. ®3.iei 

1.1 

7.. 7.. 8. 583. 8. 583 

148.8.18.8 

372.8.7.8 

8.8.8.9.148.8.115.8 

323.8.228.8 
8.34,7.8,1,2 
5.563,5.563 

1 .47,1.47 

19.75.8.125.54.8 

8.385.8.385.55.8 

8.167.18.8.18.8 
3.833,7.58 
55.8,11.18,8.819 
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A. 2 A Listing  of  SDTA2FJ  Input  File  for  Program  DBJACKS 


1.41.57.42.1 

2.42.57.58.1 

3.42.58.43.1 

4.43.58.44.1 

5.44.58.59.1 

6.44.59.45.1 

7.52.67.53.1 

8.52.66.67.1 

9.51.66.52.1 

16.56.66.51.1 

11.56.65.66.1 

12.49.65.56.1 

13.67.68.69.1 

14.67.69.53.1 

15.53.69.54.1 

16.54.69.76.1 

17.54.76.55.1 

18.55.76.56.1 

19.56.76.71.1 

26.56.71.41.1 

21.41.71.57.1 

22.57.71.72.1 

23.59.66.61.1 

24.45.59.61.1 

25.45.61.46.1 

26.46.61.62.1 

27.46.62.47.1 

28.47.62.48.1 

29.48.62.63.1 

36.48.63.49.1 

31.49.63.65.1 

32.63.64.65.1 

33.57.72.73.6 

34.57.73.74.6 

35.57.74.58.6 

36.58.74.75.6 

37.58.75.59.6 

38.59.75.76.6 

36.59.76.66.6 

46.67.84.68.6 

41.67.83.84.6 

42.66.83.67.6 

43.66.82.83.6 

44.65.82.66.6 

45.65.81.82.6 

46.64.81.65.6 

47.66.76.77.6 

48.66.77.61.6 

49.61.77.78.6 

56.61.78.62.6 

51.62.78.79.6 

52.62.79.63.6 

53.63.79.86.6 

54.63.86.64.6 

55.64.86.81.6 

56.17.74.73.6 

57.17.18.74.6 

58.18.75.74.6 

59.18.19.75.6 

66.19.26.75.6 

61.26.76.75.6 

62.26.21.76.6 

63.21.77.76.6 

64.21.22.77.6 

65.22.23.77.6 

66.31.84.83.6 

67.31.83.36.6 

68.36.83.82.6 

69.29.36.82.6 

76.28.29.82.6 

71.28.82.81.6 

72.27.28.81.6 

73.27.81.86.6 
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74.26.27.80.6 

75.25.26.80.6 

76.23.78.77.6 

77.23.24.78.6 

78.24.79.78.6 

79.24.25.79.6 

80.25.80.79.6 

81 .32.84.31 .6 

82.32.33.84.6 

83.33.34.84.6 

84.34.68.84.6 

85.34.69.68.6 

86.34.35.69.6 

87.35.36.69.6 

88.36.70.69.6 

89.36.71 .70.6 

90.36.37.71 .6 

91 .37.38.71.6 

92.38.72.71.6 

93.38.73.72.6 

94.38.39.73.6 

95.39.40.73.6 

96.40.17.73.6 

97.1 .85.86.2 

98.1 .86.2.2 

99.2.86.3.2 

100.3.86.87.2 

101 .3.87.88.2 

102.4.3.88.2 

103.5.4.88.2 

104.5.88.89.2 

105.5.89.90.2 

106.5.90.6.2 

107.6.90.7.2 

108.7.90.91 .2 

109.7.91 .92.2 

110.7.92.8.2 

111.8.92.1.2 

112.1 .92.85.2 

113.9.93.94.2 

114.10.9.94.2 

115.10.94.11.2 

116.11 .94.95.2 

117.11 .95.96.2 

118.11 .96.12.2 

119.12.96.13.2 

120.13.96.97.2 

121.13.97.98.2 

122.13.98.14.2 

123.14.98.15.2 

124.15.98.99.2 

125.15.99.100.2 

126.16.15.100.2 

127.16.100.9.2 

128.9.100.93.2 

129.85.102.86.3 

130.86.102.103.3 

131 .86.103.87.3 

132.87.103.117.3 

133.87.117.88.3 

134.103.104.117.3 

135.104.105.117.3 

136.89.117.105.3 

137.88.117.89.3 

138.89.105.90.3 

139.90.105.106.3 
140,90,106.91  .3 

141 .91 . 119.92.3 

142.92.119.85.3 

143.93.101 .94.3 

144.94.101.102.3 

145.95.94.102.3 

146.95.102.85.3 
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147.95.85.119.3 

148.96.95.119.3 

149.96.119.97.3 

150.97.119.91 .3 

151 .97.91 .106.3 

152.97.106.98.3 

153.98.106.107.3 

154.98.107.99.3 

155.99.107.118.3 

156.99.118.100.3 

157.107.108.118.3 

158.108.101 .118.3 

159.100.118.93.3 

160.93.118.101 .3 

161 .101 .109.110.4 

162.101.110.102.4 

163.102.110.103.4 

164.110.111.103.4 

165.103.111 .112.4 

166.103.112.104.4 

167.104.112.105.4 

168.112.113.105.4 

169.105.113.114.4 

170.105.114.106.4 

171 .106.114.107.4 

172.107.114.115.4 

173.107.115.116.4 

174.107.116.108.4 

175.108.116.101 .4 

176.101 .116.109.4 

177.109.128.129.5 

178.109.129.110.5 

179.110.129.130.5 

180.110.130.131 .5 

181.111.110.131 .5 

182.111.131.120.5 

183.111.120.112.5 

184.112.120.121 .5 

185.112.121 .122.5 

186.112.122.113.5 

187.113.122.123.5 

188.113.123.114.5 

189.114.123.124.5 

190.114.124.125.5 

191.114.125.115.5 

192.115.125.126.5 

193.115.126.116.5 

194.116.126.127.5 

195.116.127.128.5 

196.116.128.109.5 

197.49.129.128.5 

198.49.50.129.5 

199.50.51.129.5 

200.51.130.129.5 

201 .51 .131 .130.5 

202.51 .52.131 .5 

203.52.53.131 .5 

204.53.120.131 .5 

205.53.54.120.5 

206.54.55.120.5 

207.55.121.120.5 
208,55.122,121  ,5 

209.55.56.122.5 

210.56.41.122.5 

211 .41.123.122.5 

212.41 .42.123.5 

213.42.43.123.5 

214.43.124.123.5 

215.43.125.124.5 

216.43.44.125.5 

217.44.45.125.5 

218.45.126.125.5 

219.45.46.126.5 
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220.46.47.126.5 

221 .47.127.126.5 

222.47.128.127.5 

223.47.48.128.5 

224.48.49.128.5 
8 

82.0. 50.0.0.0.0.46.0.46 

83.0. 50.0.0.0.0.46.0.46 

86.0. 50.0.0.0.0.46.0.46 

87.0. 50.0.0.0.0.46.0.46 

90.0. 50.0.0.0.0.46.0.46 

91 .0. 50.0.0.0.0.46.0.46 

94.0. 50.0.0.0.0.46.0.46 

95.0. 50.0.0.0.0.46.0.46 


0,46.0 

0,46.0 

0,46.0 

0,46.0 

0,46.0 

0,46.0 

0,46.0 

0,46.0 
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A. 3 Output  File  SOUTFJ  from  Program  DBJACKS 


HOT  WATER  SYSTEM(P1PES  IN  ONE  CASING)AT  FORT  JACKSON,  1/7/87 
TP1  TP2  KI  KC  01  02 


323.00 

220.00 

.34 

7.00 

5.56 

5.56 

THI1 

THI2 

OIAC 

THKC 

OEPC 

SI 

S2 

TC 

1.47 

1.47 

19.75 

.13 

54.00 

.39 

.39 

55.00 

E 

m 

HY 

MONTH 

.17 

10.00 

10.00 

1 

W 

H 

0 

F 

A 

8 

WM 

HY 

5.00 

8.83 

.58 

.58 

3.83 

7.50 

10.00 

10.00 

XC1 

YC1 

XC2 

YC2 

81 

82 

2.500 

4.115 

2.500 

4.665 

4.115 

4.885 

on 

012 

S1 

S2  THI1 

KII 

KIG 

TP1 

TP2 

5.56 

5.56 

.39 

39  1.47 

.34 

7.00 

323. 

220. 

01 

02 

OT 

KP 

76.12 

36.62 

112.94 

.522 

PRANTL-  .7806  ORASOF-  .6752E48B 

CONKA-  8.8033  WADKA^  1.4159  (BTU-I»0+-rr**2-DEF  F) 


2.66 

2.73 

2.66 

2.50 

2.34 

2.27 

2.34 

2.50 

2.66 

2.50 

2.34 

2.27 

2.34 

2.50 

-10.00 

-16.00 

-16.00 

-10.00 

.00 

2.50 

5.00 

15.00 

15.00 

15.00 

15.00 

15.00 

10.00 

5.00 

4.42 

2.50 

.58 

.00 

.58 

.58 

.56 

.58 

.58 

1.54 

2.50 

3.46 

4.42 

4.42 

4.42 

3.46 

2.50 

1.54 

.00 

.06 

.58 

2.50 

4.42 

5.00 

5.00 

5.00 

5.00 

5.00 

.58 

.00 

-5.06 

-5.00 

-5.00 

-5.00 

.00 

1.67 

10.00 

10.00 

10.06 

16.00 

2.75 

2.85 

2.75 

2.50 

2.25 

2.56 

2.75 

2.85 

2.75 

2.50 

2.25 

2.15 

3.06 

3.32 

3.08 

2.50 

1.92 

1.68 

1.02 

2.50 

3.09 

2.50 

1.91 

1.67 

1.01 

2.56 

2.50 

2.50 

2.50 

1.54 

1.12 

1.12 

1.12 

1.54 

2.56 

3.46 

3.87 

Y(M).ly^1.NN 
4.28  4.12 

3.95 

3.88 

3.05 

4.12 

4.28 

4.35 

4.72 

4.65 

4.72 

4.69 

5.05 

5.12 

.25 

2.63 

13.83 

18.83 

18.83 

18.63 

18.63 

16.63 

13.63 

8.83 

.25 

.00 

.00 

.00 

.06 

.00 

.06 

.00 

.75 

2.63 

4.50 

6.36 

8.25 

6.25 

6.25 

6.25 

4.50 

2.63 

.75 

.75 

.75 

.75 

.75 

4.50 

6.83 

8.63 

8.83 

8.83 

8.25 

4.50 

.75 

.17 

.17 

.17 

.17 

.75 

4.56 

6.83 

13.83 

13.83 

8.83 

4.50 

.75 

.17 

4.37 

4.12 

3.86 

3.76 

4.37 

4.47 

5.14 

4.89 

4.63 

4.53 

4.63 

4.89 

5.08 

4.50 

3.92 

3.68 

3.02 

4.50 

5.08 

5.32 

3.91 

3.67 

3.91 

4.56 

5.09 

5.33 

.84 

.04 

2.21 

2.21 

2.63 

4.50 

6.38 

6.79 

6.79 

6.79 

2.86 

-10.00 

15.00 

-5.00 

4.42 

.00 

4.42 

3.33 

2.25 

2.25 

3.09 

2.50 

3.87 


5.05 

4.50 

4.50 

.00 

8.25 

6.25 
.17 

13.83 

3.86 

5.14 

5.09 

4.50 

6.38 


2.63 
00 

I Aft 

.00000E400 

.00000E400 

.00000E400 

.00000E400 

.00000E400 

.00000E400 

.44080E402 

.00000E400 

.6000eE400 

.000eeE400 

.0000eE400 

.0000eE400 

.0e000E406 

.000e0E400 

.445e7E402 

.44507E402 

.16136E402 

.16136E402 

.12353E402 

.00000E400 

.00000E400 

.00000E400 

.00000E40e 


ARRAY 

• WWUL  I vU 

.0000eE400 
.00000E400 

AAAAAr^.AA 

• WvvvCTW 

.0000eE40e 
.00000E40e 
.32296E403 
.26746E402 
.00006E400 
.00000E400 
.60e0eE400 
.00e0eE40e 
.00e00E400 
.0000eE40e 
.0000eE40e 
.45673E402 
.22941E402 
.1B136E402 
.18136E402 
.12353E402 
.00000E40e 
.0000eE40e 
.00e00E400 
.00000E400 


• WWOCTW 

• WI9W&T  Vv 

• WWvLI  W 

.00000E400 

.11500E403 
.64208E402 
.00000E400 
.00000E400 
.00000E400 
.00000E400 
.0000eE400 
.00000E400 
.17716E401 
. 10771 E402 
.11762E402 
.18136E402 
.12353E402 
.12353E402 
.00000E400 
.00000E400 
.0000eE400 
.00000E400 


• WVWL^W 
AAAAAr^AA 

.00000E400 

AAAAAP.XAA 

.00900E400 

.64208E402 

.11500E403 

AAAAAC..i.AA 
« WvOwb  V Uv 

.00000E400 

.00000E400 

.00000E400 

.00000E40e 

AAAAATu-AA 

• WVWbV  W 

.11762E402 

.10771E402 

.17716E+01 

.18136E402 

.12353E402 

.123S3E402 

.00000E400 

.00000E400 

AAAAAP.aAA 

AAAAAP^^IA 

• WWOL I UV 


• WVUUbTUV 

• VWWL  Wv 

• WVWCiTW 

AAAAAP.A.flkA 

• VWWbTVv 

AAAAAP.i.AA 

• VWWLf  W 

• WUUV&t  VV 

.26748E402 

.32296E403 

• WVWIC.^W 

AAAAAP.A.AA 

• WOWLTW 

AAAAAP.J.AA 

• WWv&i  W 

• WWv&TVV 

AAAAAP^A-AA 

• WWVLi  W 

.22941 E402 
.45673E402 
.1B136E402 
.1B136E402 
.12353E402 
.12353E402 

AAAAAPuAA 

.00000E400 

• WwVUCi^W 

• VWWbTVw 
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2.73 

10.00 

15.00 

10.00 

4.42 

.00 

2.50 

5.00 

2.15 

2.50 

3.33 

3.46 

3.87 


4.89 

6.83 

2.63 

.00 

6.38 

8.63 
.17 

13.83 

4.12 

5.24 

4.50 

2.21 

4.50 


.00000E+e0  .00000E+00  .00000E+00  .00000E+00  .00000E+00 
.00000E+00  .00000E-)-00  .00000E+00  .00000E+00  .00000E+00 
. 00000E+00 

PRANTL-  .7021  GRASOF-  .4822E+08 

CONKA-  7.9538  RADKA=  1.4421  (BTU-IN/H-FT**2-DEF  F) 
KASP-  9.3960  (BTU-IN./H-FT**2-DEG  F) 

AVERAGE  VALUES  OF  PIPE  INSULATION  THERMAL  CONDUCTIVITY  : 

KI1  .340  KI2  -=  .340  BTU-IN/H-FT**2-DEG  F 


AVERAGE  TEMPERATURE  DROPS  ACROSS  INSULATION 
T1-  187.05  T2>=  87.73  DEG  F 


HEAT  LOSSES  FROM  UNDERGROUND  PIPES 


01* 


02- 


36.81 


QT- 


115.30  BTU/H-FT 


78.49 
00  ARRAY 

.0e000E+00  .00000E+00  .00000E-t-00  .00000E+00  .00000E+00 

.0e000E+00  .00000E-I-00  .00000E-t-00  .00000E+00  .00000E+e0 

.00000E-f00  .00000E+00  .00000E-f00  .00000E+00  .00000E+e0 

.00000E+00  .00000E+00  .00000E+00  .00000E+00  .00000E+00 

.00000E+00  .00000E+00  .00000E-)-00  .00000E-t-00  .00000E+00 

.00000E+00  .00000E+00  .00000E-f00  .00000E-»-00  .00000E+00 

.00000E+00  .32296E+03  .11500E+03  .64208E+02  .28748E+02 

.44080E+02  .28748E+02  . 64208E+02  .11500E+03  .32296E+03 

.00000E+00  .00000E+00  .00000E+00  .00000E+00  .00000E+00 

.00000E+00  .00000E+00  .00000E+00  .00000E400  .00000E-t-00 

.00000E+00  .00000E+00  .00000E400  .00000E+00  .00000E+00 

.00000E+00  .00000E-I-00  .00000E-t-00  .00000E+00  .00000E+00 

.00000E+00  .00000E+00  .00000E+00  .00000E+00  .00000E+00 

.00000E+00  .00000E+00  .00000E+00  .00000E+00  .00000E+00 

.00000E+00  .00000E+00  .17716E+01  .11762E+02  . 22941 E+02 

.44507E+02  .45673E+02  .10771 E+02  .10771 E+02  .45673E+02 

.44507E+02  .22941 E+02  .11762E+02  .17716E+01  .16226E+02 

.17420E+02  .16019E+02  .17350E+02  .16019E+02  .17420E+02 

.16226E+02  .17598E+02  .10343E+02  .10751 E+02  .10484E+02 

.10873E+02  .10484E+02  .10751 E+02  .10343E+02  .10711 E+02 

.00000E+00  .00000E+00  .00000E+00  .00000E+00  .00000E+00 

.00000E+00  .00000E+00  .00000E+00  .00000E+00  .00000E+00 

.00000E+00  .00000E+00  .00000E+00  .00000E+00  .00000E+00 

.00000E+00  .00000E+00  .00000E+00  .00000E+00  .00000E+00 

.00000E+00  .00000E+00  .00000E+00  .00000E+00  .00000E+00 

.00000E+00  .00000E+00  .00000E+00  .00000E+00  .00000E+00 

.00000E+00 

PRANTL-  .7028  GRASOF-  .5110E+08 

CONKA«  8.0130  RADKA«=  1.4244  (BTU-IN/H-FT**2-DEF  F) 
KASP-  9.4375  (BTU-IN ./H-FT**2-DEG  F) 


AVERAGE  VALUES  OF  PIPE  INSULATION  THERMAL  CONDUCTIVITY  : 

KI1  ■ .318  KI2  - .293  BTU-IN/H-FT**2-DEG  F 


AVERAGE  TEMPERATURE  DROPS  ACROSS  INSULATION  : 
T1-  191.57  T2-  92.47  DEG  F 


HEAT  LOSSES  FROM  UNDERGROUND  PIPES 


01-  75.17  02- 

00  ARRAY 
.00000E+00  .00000E+00 


.00000E+00 
.00000E+00 
. 00000E+00 
.00000E+00 
.00000E+00 
. 00000E+00 
. 44080E+02 
. 00000 E+00 
.00000E+00 
.00000E+00 
.00000E+00 
.00000E+00 
.00000E+00 
.00000E+00 
.44507E+02 
.44507E+02 


.00000E+00 
. 00000E+00 
.00000E+00 
.00000E+00 
. 00000E+00 
.32296E+03 
. 28748E+02 
. 00000E+00 
. 00000E+00 
.00000E+00 
. 00000E+00 
. 00000E+00 
. 00000E+00 
. 00000E+00 
.45673E+02 
.22941 E+02 


33.42  OT- 

. 00000E+00 
.00000E+00 
. 00000E+00 
.00000E+00 
.00000 E+00 
.00000E+00 
.11500E+03 
. 64208E+02 
.00000E+00 
.00000 E+00 
.00000 E+00 
.00000 E+00 
.00000 E+00 
.00000 E+00 
. 17716E+01 
.10771 E+02 
. 11762E+02 


108.59 

. 00000E+00 
.00000E+00 
. 00000E+00 
.00000E+00 
. 00000E+00 
.00000E+00 
. 64208E+02 
. 11500E+03 
.00000 E+00 
.00000E+00 
.00000E+00 
.00000E+00 
. 00000E+00 
. 00000E+00 
.11762E+02 
.10771 E+02 
. 17716E+01 


BTU/H-FT 

.00000 E+00 
.00000E+00 
.00000 E+00 
.00000 E+00 
. 00000E+00 
. 00000E+00 
. 28748E+02 
.32296E+03 
.00000E+00 
.00000 E+00 
.00000E+00 
.00000 E+00 
.00000 E+00 
.00000 E+00 
.22941 E+02 
. 45673E+02 
.16156E+02 
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.17371E+e2 
.16156E+02 
.10848E-f02 
. 00000 E-hCe 
.00000E+00 
. 00000 E+00 
. 00000E+00 
. 00000 E-t-00 
. 00000E+00 
.00000E+00 
PRANTb 
CONKA« 
KASP- 


.15945E+02 
. 17554E+02 
, 1 0440E+02 
.00000 E+00 
.00000E+00 
. 00000 E+0e 
.00000 E+00 
.00000 E+00 
.00000 E+00 


.17298E+02 
.10296E+02 
. 10721 E+02 
. 00000E+00 
.00000 E+00 
. 00000E+00 
.00000 E+00 
.00000 E+00 
.00000 E+00 


.15945E+02 
.10721 E+02 
.10296E+02 
. 00000E+00 
.00000 E+00 
.00000E+00 
.00000 E+00 
.00000 E+00 
.00000 E+00 


.17371 E+02 
. 1 0440E+02 
. 1 0680E+02 
. 00000E+00 
. 00000E+00 
.00000 E+00 
. 00000E+00 
. 00000E+00 
. 00000E+00 


.7028 
8.0017 
9 . 4257 


GRASOF-  .5086E+08 

RADKA=  1.4239  (BTU-IN/H-FT**2-0EF  F) 
(BTU-IN./H-FT**2-DEG  F) 


AVERAGE  VALUES  OF  PIPE  INSULATION  THERMAL  CONDUCTIVITY  : 

KI1  - .317  KI2  ■=  .292  BTU-IN/H-FT**2-DEG  F 


AVERAGE  TEMPERATURE  DROPS  ACROSS  INSULATION  : 
T1-  191.80  T2-  92.67  DEG  F 


HEAT  LOSSES  FROM  UNDERGROUND  PIPES  : 

01-  75.00  02-  33.38  OT-  108.38  BTU/H-FT 

M I J K MAT . C 


1 

41 

57 

42 

1 

.5833 

2 

42 

57 

58 

1 

.5833 

3 

42 

58 

43 

1 

.5833 

4 

43 

58 

44 

1 

.5833 

5 

44 

58 

59 

1 

.5833 

6 

44 

59 

45 

1 

.5833 

7 

52 

67 

53 

1 

.5833 

8 

52 

66 

67 

1 

.5833 

9 

51 

66 

52 

1 

.5833 

10 

50 

66 

51 

1 

.5833 

11 

50 

65 

66 

1 

.5833 

12 

49 

65 

50 

1 

.5833 

13 

67 

68 

69 

1 

.5833 

14 

67 

69 

53 

1 

.5833 

15 

53 

69 

54 

1 

.5833 

16 

54 

69 

70 

1 

.5833 

17 

54 

70 

55 

1 

.5833 

18 

55 

70 

56 

1 

.5833 

19 

56 

70 

71 

1 

.5833 

20 

56 

71 

41 

1 

.5833 

21 

41 

71 

57 

1 

.5833 

22 

57 

71 

72 

1 

.5833 

23 

59 

60 

61 

1 

.5833 

24 

45 

59 

61 

1 

.5833 

25 

45 

61 

46 

1 

.5833 

26 

46 

61 

62 

1 

.5833 

27 

46 

62 

47 

1 

.5833 

28 

47 

62 

48 

1 

.5833 

29 

48 

62 

63 

1 

.5833 

30 

48 

63 

49 

1 

.5833 

31 

49 

63 

65 

1 

.5833 

32 

63 

64 

65 

1 

.5833 

33 

57 

72 

73 

6 

.5833 

34 

57 

73 

74 

6 

.5833 

35 

57 

74 

58 

6 

.5833 

36 

58 

74 

75 

6 

.5833 

37 

58 

75 

59 

6 

.5833 

38 

59 

75 

76 

6 

.5833 

39 

59 

76 

60 

6 

.5833 

40 

67 

84 

68 

6 

.5833 

41 

67 

83 

84 

6 

.5833 

42 

66 

83 

67 

6 

.5833 

43 

66 

82 

83 

6 

.5833 

44 

65 

82 

66 

6 

.5833 

45 

65 

81 

82 

6 

.5833 

46 

64 

81 

65 

6 

.5833 

47 

60 

76 

77 

6 

.5833 

48 

60 

77 

61 

6 

.5833 

49 

61 

77 

78 

6 

.5833 

50 

61 

78 

62 

6 

.5833 

47 


51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

110 

1 1 1 

112 

113 

114 

115 

116 

117 

118 

119 

120 

121 

122 

123 


62 

78 

79 

6 

.5833 

62 

79 

63 

6 

.5833 

63 

79 

80 

6 

.5833 

63 

80 

64 

6 

.5833 

64 

80 

81 

6 

.5833 

17 

74 

73 

6 

.5833 

17 

18 

74 

6 

.5833 

18 

75 

74 

6 

.5833 

18 

19 

75 

6 

.5833 

19 

20 

75 

6 

.5833 

20 

76 

75 

6 

.5833 

20 

21 

76 

6 

.5833 

21 

77 

76 

6 

.5833 

21 

22 

77 

6 

.5833 

22 

23 

77 

6 

.5833 

31 

84 

83 

6 

.5833 

31 

83 

30 

6 

.5833 

30 

83 

82 

6 

.5833 

29 

30 

82 

6 

.5833 

28 

29 

82 

6 

.5833 

28 

82 

81 

6 

.5833 

27 

28 

81 

6 

.5833 

27 

81 

80 

6 

.5833 

26 

27 

80 

6 

.5833 

25 

26 

80 

6 

.5833 

23 

78 

77 

6 

.5833 

23 

24 

78 

6 

.5833 

24 

79 

78 

6 

.5833 

24 

25 

79 

6 

.5833 

25 

80 

79 

6 

.5833 

32 

84 

31 

6 

.5833 

32 

33 

84 

6 

.5833 

33 

34 

84 

6 

.5833 

34 

68 

84 

6 

.5833 

34 

69 

68 

6 

.5833 

34 

35 

69 

6 

.5833 

35 

36 

69 

6 

.5833 

36 

70 

69 

6 

.5833 

36 

71 

70 

6 

.5833 

36 

37 

71 

6 

.5833 

37 

38 

71 

6 

.5833 

38 

72 

71 

6 

.5833 

38 

73 

72 

6 

.5833 

38 

39 

73 

6 

.5833 

39 

40 

73 

6 

.5833 

40 

17 

73 

6 

.5833 

1 

85 

86 

2 

.0251 

1 

86 

2 

2 

.0277 

2 

86 

3 

2 

.0277 

3 

86 

87 

2 

.0249 

3 

87 

88 

2 

.0249 

4 

3 

88 

2 

.0276 

5 

4 

88 

2 

.0276 

5 

88 

89 

2 

.0249 

5 

89 

90 

2 

.0249 

5 

90 

6 

2 

.0277 

6 

90 

7 

2 

.0277 

7 

90 

91 

2 

.0251 

7 

91 

92 

2 

.0254 

7 

92 

8 

2 

.0280 

8 

92 

1 

2 

.0280 

1 

92 

85 

2 

.0254 

9 

93 

94 

2 

.0237 

10 

9 

94 

2 

.0248 

10 

94 

11 

2 

.0248 

11 

94 

95 

2 

.0238 

11 

95 

96 

2 

.0241 

11 

96 

12 

2 

.0251 

12 

96 

13 

2 

.0251 

13 

96 

97 

2 

.0241 

13 

97 

98 

2 

.0238 

13 

98 

14 

2 

.0248 

14 

98 

15 

2 

.0248 
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15 

15 

16 

16 

9 

85 

86 

86 

87 

87 

103 

104 

89 

88 

89 

90 

90 

91 

92 

93 

94 

95 

95 

95 

96 

96 

97 

97 

97 

98 

98 

99 

99 

107 

108 

100 

93 

101 

101 

102 

110 

103 

103 

104 

112 

105 

105 

106 

107 

107 

107 

108 

101 

109 

109 

110 

110 

111 

111 

111 

112 

112 

112 

113 

113 

114 

114 

114 

115 

115 

116 

116 

116 


98 

99 

2 

.0237 

99 

100 

2 

.0236 

15 

100 

2 

.0247 

100 

9 

2 

.0247 

100 

93 

2 

.0236 

102 

86 

3 

.7865 

102 

103 

3 

.7865 

103 

87 

3 

.7865 

103 

117 

3 

.7865 

117 

88 

3 

.7865 

104 

117 

3 

.7865 

105 

117 

3 

.7865 

117 

105 

3 

.7865 

117 

89 

3 

.7865 

105 

90 

3 

.7865 

105 

106 

3 

.7865 

106 

91 

3 

.7865 

119 

92 

3 

.7865 

119 

85 

3 

.7865 

101 

94 

3 

.7865 

101 

102 

3 

.7865 

94 

102 

3 

.7865 

102 

85 

3 

.7865 

85 

119 

3 

.7865 

95 

119 

3 

.7865 

119 

97 

3 

.7865 

119 

91 

3 

.7865 

91 

106 

3 

.7865 

106 

98 

3 

.7865 

106 

107 

3 

.7865 

107 

99 

3 

.7865 

107 

118 

3 

.7865 

118 

100 

3 

.7865 

108 

118 

3 

.7865 

101 

118 

3 

.7865 

118 

93 

3 

.7865 

118 

101 

3 

.7865 

109 

110 

4 

31 .0000 

110 

102 

4 

31 .0000 

110 

103 

4 

31 .0000 

111 

103 

4 

31 .0000 

111 

112 

4 

31 .0000 

112 

104 

4 

31 .0000 

112 

105 

4 

31 .0000 

113 

105 

4 

31 .0000 

113 

114 

4 

31 .0000 

114 

106 

4 

31 .0000 

114 

107 

4 

31 .0000 

114 

115 

4 

31 .0000 

115 

116 

4 

31 .0000 

116 

108 

4 

31 .0000 

116 

101 

4 

31 .0000 

116 

109 

4 

31 .0000 

128 

129 

5 

.5833 

129 

110 

5 

.5833 

129 

130 

5 

.5833 

130 

131 

5 

.5833 

110 

131 

5 

.5833 

131 

120 

5 

.5833 

120 

112 

5 

.5833 

120 

121 

5 

.5833 

121 

122 

5 

.5833 

122 

113 

5 

.5833 

122 

123 

5 

.5833 

123 

114 

5 

.5833 

123 

124 

5 

.5833 

124 

125 

5 

.5833 

125 

115 

5 

.5833 

125 

126 

5 

.5833 

126 

116 

5 

.5833 

126 

127 

5 

.5833 

127 

128 

5 

.5833 

128 

109 

5 

.5833 

49 


197 

49 

129 

128 

5 

.5833 

198 

49 

50 

129 

5 

.5833 

199 

50 

51 

129 

5 

.5833 

200 

51 

130 

129 

5 

.5833 

201 

51 

131 

130 

5 

.5833 

202 

51 

52 

131 

5 

.5833 

203 

52 

53 

131 

5 

.5833 

204 

53 

120 

131 

5 

.5833 

205 

53 

54 

120 

5 

.5833 

206 

54 

55 

120 

5 

.5833 

207 

55 

121 

120 

5 

.5833 

208 

55 

122 

121 

5 

.5833 

209 

55 

56 

122 

5 

.5833 

210 

56 

41 

122 

5 

.5833 

211 

41 

123 

122 

5 

.5833 

212 

41 

42 

123 

5 

.5833 

213 

42 

43 

123 

5 

.5833 

214 

43 

124 

123 

5 

.5833 

215 

43 

125 

124 

5 

.5833 

216 

43 

44 

125 

5 

.5833 

217 

44 

45 

125 

5 

.5833 

218 

45 

126 

125 

5 

.5833 

219 

45 

46 

126 

5 

.5833 

220 

46 

47 

126 

5 

.5833 

221 

47 

127 

126 

5 

.5833 

222 

47 

128 

127 

5 

.5833 

223 

47 

48 

128 

5 

.5833 

224 

48 

49 

128 

5 

.5833 

TEMPERATURE 

ARRAY  : 

T(I). 

1=1 .NN 

.32300E+03 
.32300E+03 
.22000E+03 
.22000E+03 
. 54680E+02 
. 55395E+02 
.45533E+02 
.58719E+02 
.64637E+02 
. 80951 E+02 
.95188E+02 
.66091 E+02 
.76850E+02 
.86835E+02 
.59 157 E+02 
.62056E+02 
.62056E+02 
. 12995E+03 
. 13652E+03 
. 1 4728E+03 
. 1 1 648E+03 
.11728E+03 
. 1 1898E+03 
.11521E+03 
.85590E+02 
.90360E+02 
.85938E+02 


. 32300E+03 
.32300E+03 
.22000E+03 
. 45533E+02 
.55395E+02 
.54680E+02 
.45423E+02 
.57478E+02 
.80731 E+02 
.81831 E+02 
.80731 E+02 
.63715E+02 
.78578 E+02 
.63715E+02 
.59022E+02 
. 62748E+02 
.60832 E+02 
.12256E+03 
. 14814E+03 
.13369E+03 
.11728E+03 
. 1 1 648E+03 
. 11717E+03 
.12109E+03 
.82806E+02 
.92821 E+02 


.32300E+03 
.32300E+03 
. 22000E+03 
. 47363E+02 
. 55395E+02 
.52456E+02 
.48811 E+02 
.5771 5 E+02 
.95188E+02 
.80951 E+02 
. 64637E+02 
. 86835E+02 
. 76850E+02 
.59022 E+02 
.49342 E+02 
.63459 E+02 
.51146E+02 
. 12338E+03 
. 11803 E+03 
. 12476E+03 
, 11898E+03 
. 11523E+03 
.11898 E+03 
.11 750 E+03 
. 85938E+02 
. 90360E+02 


. 32300E+03 
.22000E+03 
. 22000E+03 
. 49045E+02 
.55395E+02 
.49045E+02 
.5771 5 E+02 
.48811 E+02 
.87191E+02 
. 79263E+02 
.66091 E+02 
.779 50 E+02 
. 75833E+02 
.59157E+02 
.51146E+02 
. 63459E+02 
.49 342 E+02 
.12256E+03 
. 12476E+03 
.11803E+03 
.11719E+03 
.11647 E+03 
.11728E+03 
.14756E+03 
.10489 E+03 
.91 640 E+02 


. 32300 E+03 
.22000 E+03 
. 22000E+03 
.52456E+02 
.55395E+02 
.47363E+02 
. 57478E+02 
.45423 E+02 
.79263E+02 
.871 91 E+02 
. 66905E+02 
.75833 E+02 
.77950E+02 
.60490E+02 
. 60832E+02 
. 62748E+02 
.13652 E+03 
.12995E+03 
. 13369E+03 
. 1 1 836E+03 
.11898E+03 
.11728 E+03 
. 1 1 647 E+03 
.82806E+02 
.91 640 E+02 
.10489 E+03 


AVERAGE  VALUES  OF  PIPE  INSULATION  THERMAL  CONDUCTIVITY  : 

KI1  -=  .317  KI2  - .292  BTU-IN/H-FT**2-DEG  F 


AVERAGE  TEMPERATURE  DROPS  ACROSS  INSULATION  ; 

T1-=  191  .80  T2-  92.67  DEG  F 

HEAT  LOSSES  FROM  UNDERGROUND  PIPES  : 

Q1«=  75.00  Q2-  33.38  OT-  108.38  BTU/H-FT 
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APPENDIX  B.  A Listing  of  Computer  Programs 


PROGRAM  DIRECTI 

C THIS  IS  A MAIN  PROGRAM  FOR  HEAT  LOSS  ANALYSIS  OF  DIRECTLY  BURIED 
C CONDUIT  UNDERGROUND  HEAT  DISTRIBUTION  SYSTEMS  WITH  BOTH  INSULATED 
C PIPES  IN  THE  SAME  CASING.  BASED  ON  THE  FINITE  ELEMENT  METHOD  USING 
C THREE  - NODE  LINEAR  TRIANGULAR  ELEMENTS. 

C SUBROUTINES  CALLED:  PIPED. TOO. SOILK, INSULK. TCXX.SOLVLE.PIPEHL.TWOPIP, 
C EOUIKO. 

C INPUT  DATA  FILES:  SOATA1  AND  SOATA2 

C OUTPUT  FILE:  SOUTPUT 

C X(n:  THE  X-COORDINATE  OF  NODAL  POINT  I. IN  FT 

C Y(I):  THE  Y-COORDINATE  OF  NODAL  POINT  I. IN  FT 

C (NOOE(M.I).I-1.3):  THREE  NODAL  POINTS  OF  ELEMENT  M 
C M ELEMENT  INDEX 

C NE  TOTAL  NUMBER  OF  ELEMENTS 

C NN  TOTAL  NUMBER  OF  NODAL  POINTS 

C MZ  TOTAL  NUMBER  OF  KNOWN  NODAL  TEMPERATURES 

C C THERMAL  CONDXTIVITY.  BTU-IKA«/FT**2/DEG  F 

C L THICKNESS  OF  THE  ELEMENT.  H 

C T(I):  THE  TEMPERATURE  OF  NODAL  POINT  I.  IN  DEC  F 
REAL  L.KK.KI .KG.KIX1 ,KIX2.KTCT.KASP.KCAS.KBF 
CHARACTER*4  TITLE(15) 

DIMENSION  Q(15e).T(15e).X(15e).Y(15e).KK(15e.15e) 

DIMENSION  AS(25e) .B2IZ(25e) .B3IZ(25e) .B2JZ(25e) .B2KZ(25e) . 
k B3JZ(25e).B3KZ(25e) 

DIMENSION  CC(25e) .TCX(12.6) .QQ(15«) .NOOE(25e.3) .MAT(25e) 
DIMENSION  HIJ(25e).HJK(25e).HKI(25e).TIJ(25e).TJK(25e). 
ft  TKI (250) .HHIJ(25e) .HHJK(25e) .HHKI (250) . IXC6(250) 

DIMENSION  CK(150. 150), 00058). XT(150).INDX(150).W(150) 
C0Mi40N/PP/TP1  .TP2.KI  .KC.D1  .D2.TH1  .TH2.DP1  .DP2.S1  .S2.TC. 
ft  WW.HY. MONTH 

COIftlON  /EKC/D1P.02P.DIC.THK1.THK2 
COMilON  /ST/AO.BO.DIFF 
PI-4.«ATAN(1.) 

OPEN  (B.FILE-’SOATAr) 

OPEN  (7.FILE-'S0UTPUT',STATUS-*NEW'.F0Rfc^*F0RMATTED') 

OPEN  (9,FILE-’S0ATA2’) 

C READ  IN  THE  TITLE  OF  THE  PROBLEM  TO  BE  ANALYZED 
READ  (8.2.ERR-20e0)  TITLE 

2 FORMAT (15A4) 

WRITE  (7.3)  TITLE 

3 FORMAT (1X.15A4) 

C READ  TOTAL  NUMBER  OF  NODAL  POINTS.  TOTAL  NUMBER  OF  TRIANGULAR 
C ELEMENTS.  TOTAL  NUifiER  OF  KNOWN  NODAL  TEkPERATURES . AND  THE 
C FIRST  ELEMENT  IM>EX  OF  PIPE  INSULATION 
READ  (6.*)  NN.NE.MZ.MINS 
C SET  THE  UNIT  NUkBER  OF  THE  PRINTER 
M07 

C READ  MONTH  OF  INTEREST  AND  THE  INDEX  FOR  FINITE  ELEMENT  GRID  DATA 
C TO  BE  PRINTED  OUT  : ICALB  - 1 PRINT  OUT  NODAL  COORDINATES 
C • B NO  PRINT  OUT 

READ  (8.*)  MONTH. ICALB 

C READ  THE  THERMAL  CONDUCTIVITY  (IN  BTU-IN./H-FT*»2  - DEC  F). 

C THICKNESS  (IN  INCHES)  OF  THE  SIDE.  THE  DEPTH  (IN  FT.)  OF  EARTH 
C COVER  (IN  FT.).  AND  THE  THICKNESS  (IN  FT.)  OF  BOTTOM  BED  OF 
C THE  INNER  EARTH  REGION. 

READ  (8.*)  KTCT.TRTK.D.F 

C READ  IN  THE  ESTIMATED  AVERAGE  TEkPERATURE  OF  AIR  INSIDE  THE 
C AIRSPACE  BETWEEN  INSULATED  PIPES  AND  OUTER  CASING.  IN  DEG  F. 

C AND  THE  TEMPERATURE  DIFFERENCE  BETWEEN  THE  INSULATED  PIPE 
C SURFACE  TEkPERATURE  AND  THE  INNER  SURFACE  TEkPERATURE  OF  THE 
C CASING.  IN  DEG  F. 

READ  (8.«)  TAS.TDEL 

C READ  IN  THE  THERMAL  CONDUCTIVITY  (IN  BTU-IN/H-FT**2-0EC  F)  OF 
C THE  PIPE  CASING  AND  THE  CONDUCTIVITY  OF  POURED-IN  INSULATION 
C MATERIAL  OR  BACK^^FILL  SOIL  SURROUNDING  THE  PIPES  IN  THE 
C INNERMOST  REGION. 

READ  (8.*)  KCAS.KBF 

C READING  IN  INPUT  DATA  FOR  CALCULATIONS  OF  PIPE  HEAT  LOSS  AND 
C GENERATION  OF  THE  COORDINATES  OF  NOOAL  POINTS 
CALL  PIPEO(X,Y.TRTK.D.F.INXK) 

CALL  TWOPIP(I) 

CALL  EQUIKO(TAS.TDEL.KASP) 

IF(ICALB.EQ.I)  THEN 
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ooooooooo 


WRITE(7.5) 

5 FORMAT(*  X(M)  .1^1  .NN’ ) 

WRITE(7.7)  (X(I).I-I .NN) 

7 FORMAT(10F7.2) 

WRITE(7, 10) 

10  FORMAT (■  Y(M) .M=1 .NN* ) 

WRITE(7.7)  (Y(I).I-I .NN) 

END  IF 

C CALCULATIONS  OF  UNDISTURBED  EARTH  TEMPERATURES  AT  VARIOUS  DEPTHS 
CALL  TGO(TGX.PI ,Y) 

C INITIALIZATION  OF  THE  INDEX  OF  CONVECTION  BOUNDARY  FOR  ELEMENT  N 
DO  12  N=1 .NE 

12  IXCB(N)-0 

C PERFORM  ITERATIONS  TO  ACCOUNT  FOR  THE  TEMPERATURE  EFFECTS  ON  SOIL 

C AND  INSULATION  THERMAL  CONDUCTIVITIES 
DO  24  I-=1  ,NN 

24  T(I)=TG 

DO  26  1=1 .NE 
HIJI 
HJK( 

HKI 
TUI 
TJK( 

TKK 
HHIJI 
HHJKI 
HHKI I 

26  CONTINUE 

READING  IN  THE  ELEMENT  NUMBER  AND  ITS  NODAL  POINTS  AND  THE 
MATERIAL  TYPE,  WHICH  INCLUDES 

MAT(J)  =1  SOIL  IN  INNER  EARTH  REGION 
= 2 PIPE  INSULATION 

= 3 AIR  SPACE  SURROUNDING  THE  PIPE  IN  CASING 
= 4 OUTER  CASING  OF  THE  INSULATED  PIPE 
- 5 POURED-IN  INSULATION  MATERIAL  IN  BACK-FILL 
REGION 

= 6 NATIVE  SOIL  IN  THE  OUTER  EARTH  REGION 


DO  30  1=1 ,NE 

READ(9,*)  J.(NODE(J.K).K=1 .3).MAT(J) 

IF  (MAT(J).EO.I)  CC(j)-KTCT/12. 

IF  (MAT(J).E0.2)  CC(J)-KI/12. 

IF  (MAT(J).EQ.3)  CC( J )-KASP/12. 

IF  (MAT(j).EQ.4)  CC( J )»KCAS/12 . 

IF  (MAT(J).EQ.5)  CC{ j)-KBF/12 . 

IF  (MAT(J).EQ.6)  CC(J)-KG/12. 

30  CONTINUE 

C READ  IN  TOTAL  NUMBER  OF  ELEMENTS  HAVING  BOUNDARY  SEGMENTS  SUBJECT 
C TO  CONVECTIVE  HEAT  TRANSFER 
READ  (9.*)  NECB 

C READ  IN  ELEMENT  NUMBER.  CONVECTIVE  HEAT  TRANSFER  COEFFICIENTS. 

C AND  AMBIENT  TEMPERATURES  FOR  THREE  BOUNDARY  SEGMENTS 
DO  35  1-1 .NECB 

READ  (9.*)  M.HIJ(M).HJK(M),HKI(M).TIJ(M).TJK(M).TKI(M) 
IXCB(M)-1 
35  CONTINUE 

ITER-1 

38  DO  40  1-1. NN 

DO  40  J-1 .NN 
Q(I)-0. 

KK(I . J)-0. 

QQ(i)-0. 

CK( I ,J)-0. 

DQ( I)-0. 

VV(I)-1 .0 
INDXO)“1 
40  CONTINUE 


L-1  . 

DO  180  M-1  .NE 
I-NODEfM. 1) 

J-NODE(M.2) 

K-N0DE(M.3) 

IF(MAT(M) .EQ.3)  CC(M)-KASP/12. 
C-CC(M) 
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IF  ((INXK.EQ.0).OR.(ITER.EQ.1))  GO  TO  60 
C DETERMINE  SOIL  AND  INSULATION  THERMAL  CONDUCTIVITIES  BASED  ON  THE 
C MEAN  TEMPERATURES 

T^*-(T(I)+T(J)+T(K))/3. 

IF(MAT(M) .EQ.2)  CALL  INSULK(TM,C) 

IF(MAT(m) .EQ.6)  CALL  SOI LKfTM.KG.C) 

IF(MAT(M).EQ.1)  CALL  SOI LK(TM.KG,C) 

CC(M)-C 
60  Xl-Xfl) 

XJ=X(j) 

XK=X(K) 

Yi-yn) 

YJ=Y( j) 

YK-=Y(K) 

CXX-C 

CXY-0, 

CYX-0 . 

CYY=C 

B2I=YJ-YK 

B3I-=XK-XJ 

B2J=YK-YI 

B3J=XI-XK 

B2K=YI-YJ 

B3K=X  J*"X  I 

C CALCULATE  THE  ELEMENT  AREA 

SA-0.5*(XJ*B2J+XI*B2I+XK*B2K) 

SA»ABS(SA) 

A2=SA*2. 

AS(M)-A2 

B2I«=B2I/A2 

B3I»*B3I/A2 

B2J=B2J/A2 

B3J=B3J/A2 

B2K«B2K/A2 

B3K-e3K/A2 

B2IZ(M)-B2I 

B3IZ(m)-B3I 

B2JZ(M)-B2J 

B3JZ(M)«B3J 

B2KZ(M)»B2K 

B3KZ(M)-B3K 

BIIi=SA*L*(B2I*B2I*CXX+B2I*B3I*CXY+B3I*B2I*CYX+B3I*B3I*CYY) 

B I J=SA*  L* { B2 1 *82 J •CXX+B2 1 •B3J •CXY+B3 1 *62 J •CYX+B3 1 •B3J  *CYY) 

B I K=SA*  L* ( B2 1 •B2K*CXX+B2 1 •B3K*CXY+B3 1 •B2K*CYX+B3 1 *B3K*CYY ) 

B J I-=SA*  L*  ( B2  J *82 1 *CXX+B2  J *83 1 •CXY+B3  J *82 1 -CYX+BSJ  *83 1 •CYY  ) 
BJJ=SA*L*(B2J*B2J*CXX+B2J*B3J*CXY+B3J*B2J*CYX+B3J*B3J*CYY) 
BJK-SA*L*(B2J*B2K*CXX+B2J*B3K*CXY+B3J*B2K*CYX+B3J*B3K*CYY) 
BKI«=SA*  L* (B2K*B2 1 •CXX+B2K*B3I •CXY+B3K*B2I •CYX+B3K*B3I •CYY) 
BKJ-SA*L*(B2K*B2J*CXX+B2K*B3J*CXY+B3K*B2J*CYX+B3K*B3J*CYY) 
BKK-SA*  L* (B2K*B2K*CXX+B2K*B3K*CXY+B3K*B2K*CYX+B3K*B3K*CYY) 


KK( 

;i, 

.1] 

l-KKI 

;i.i 

)+BII 

KK| 

[1. 

.J 

l-KKI 

I .J 

)+BIJ 

KK( 

I, 

.K 

l-KKI 

I .K 

)+BIK 

KK| 

l-KKI 

J.I 

)+BJI 

KK| 

-J 

l-KKI 

J.J 

i+BJ  J 

KKI 

J, 

.K 

l-KKI 

J.K 

Wbjk 

KK| 

K, 

.1 

l-KKI 

K.I 

)+BKI 

KKI 

K, 

.J 

l-KKI 

K.  J 

)+BKJ 

KKI 

K 

,K 

l-KKI 

[K.K 

)+BKK 

IFI 

;iXCB(M), 

,EQ.( 

0)  GO  TO  130 

C ADDITION  OF  CONVECTION  TERMS  TO  THE  ELEMENT  MATRIX  TO  ACCOUNT 
C FOR  CONVECTION  ON  BOUNDARY 

C READING  IN  CONVECTIVE  HEAT  TRANSFER  COEFFICIENTS  AND  AKfilENT 
C TEMPERATURES  FOR  THREE  BOUNDARY  SEGMENTS 

HHIJfM)-HIJ(M)*L*SQRT((Xn)-X(  J))**2+fYf  n-Y(J))**2)/6. 
HHJK(M)-HJKfM)*L*SQRTf (X(j)-XfK))**2+(Yf j)-YfK))**2)/6. 
HHKI(M)-HKI(M)*L*SORT((X(K)-X(I))**2+{Y(K)-Y(I))**2)/6. 

KK(I  .n-HHIJfM)»2.+HHKI(M)*2.+KK(I . I) 

KKfl .JI-HHIJfMI+KKfl 
KK(  I .K)-HHKI  (M)+KKn  .K) 

KKf J . l)-HHIJ(M)+KK( J . I ) 

KKf J,J)-HHIJ(M)*2+HHJK(M)*2.+KK(J.J) 

KK(J,K)«HHJK(M)+KK(J.K) 
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KKfK.  n-HHKUM^+KKrK,  n 
kk(k, j)-hhjk(m)+kk(k. J) 
KK(K.K)-HHJK(M)*2.+HHKI(M)*2.+KK(K,K) 

hhjk(m)-=tjk(m)*3.*hhjk(m) 

HHKI (M)-TKI (M)*3. •HHKI (m) 

130  Qf  n«=On)+HHIJrM)+HHKl(M) 

0(J)-Q( j)+hhij(m)+hhjk(m) 
o(k)«o(khhhjk(m)+hhki  (M) 


180 

CONTINUE 

185 

FORMAT  f ’ M I J K 

MAT . C ’ ) 

187 

FORMATOX.5I6.F10.4) 

C 

DETERMINE  OUTER  SURFACE  TEMPERATURES 
DO  200  1-1 .8 

T(I)-TP1 

II-I+8 

OF  UNDERGROUND  PIPES 

T(II)-TP2 

- 

200 

CONTINUE 

C 

DETERMINE  OUTER  BOUNDARY  TEMPERATURES 
CALL  TGXX(T.TGX. MONTH) 

MZ1-MZ+1 

DO  260  I-MZ1 .NN 

SUM-0. 

DO  250  J-1 .MZ 

OF  EARTH  REGION 

250 

SUM=SUM+KK(I . J)*T(J) 
OQ(I)-Q(I)-SUM 

260 

CONTINUE 

IF(ICALB.EQ.I)  THEN 

WRITE(7.280) 

280 

FORMAT ( 6X . • QQ  ARRAY  * ) 

WRITE(7.285)  (QQ( I ) . 1-1 .NN) 

285 

FORMAT  (5E12.5) 

END  IF 

C RENAMING  OF  MATRICES 
MN=NN-MZ 
DO  300  1-1 .MN 
K-M2+I 

DO  290  J-1 .MN 
KL-MZ+J 

290  CKfl . J)-KK(K.KL) 

XT(n=T(K) 

DQ(I)-QQ(K) 

300  CONTINUE 

C SOLUTION  OF  SIMULTANEOUS  EQUATIONS 
C SET  PHYSICAL  DIMENSION  OF  MATRIX  A 
NP-150 

CALL  SOLVLE(CK.MN.NP.  INDX.W.DO) 

DO  310  1-1 .MN 
KHWIZ+I 
T(K)-DQ(I) 

310  CONTINUE 

C CALCULATIONS  OF  EQUIVALENT  THERMAL  CONDUCTIVITIES  FOR  AIR 
C SPACE  INSIDE  THE  OUTER  CASING 
C 

C CALCULATE  AVERAGE  SURFACE  TEMPERATURE  OF  INSULATED  PIPES 
SU1-0.0 
SU2-0 . 0 
DO  312  1-1 .8 
L1-I+80 
L2-I+88 
SU1-SU1+T(L1) 

SU2-SU2+T(L2) 

312  CONTINUE 

TSM1-SU1/8. 

TSM2-SU2/8 . 

C DETERMINE  THE  EFFECTIVE  SURFACE  TEMPERATURE  OF  INSULATED  PIPES. 
C IN  DEG  F. 

DIP1-D1P+2.*THK1 

DIP2-D2P+2.*THK2 

TEFPS-(DIP1*TSM1+DIP2*TSM2)/(DIP1+DIP2) 

C CALCULATE  THE  INNER  SURFACE  TEMPERATURE  OF  THE  CASING 
SU3-0 . 0 
DO  314  1-1 .8 
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L3-I+96 

314  SU3-SU3+T(L3) 

TSIOSU3/8. 

C DETERMINE  THE  TEMPERATURE  DIFFERENCE  BETWEEN  THE  EFFECTIVE  SURFACE 
C TEMPERATURE  OF  INSULATED  PIPES  AND  THE  INNER  SURFACE  OF  THE  CASING. 
C IN  DEG  F 

TDEL-ABS(TEFPS-TSIC) 

TAS-(TEFPS+TSIC)/2. 

CALL  EQUIKO(TAS.TDEL.KASP) 

WRITE(7.320)  KASP 

320  FORMATf*  KASP=’ . F10.4.2X. ' (BTU-IN./H-FT**2-DEG  F) * ) 

330  FORMAT('  TEMPERATURE  ARRAY  : T(I).  1=1 .NN  ’) 

C CALCULATE  THE  MEAN  VALUES  OF  INSULATION  THERMAL  CONDUCTIVITY  FOR 
C PIPES  1 AND  2 
350  SK1 1-0. 

SKI2-0. 

DO  400  LN-1 ,16 
LM-MINS+LN-1 
LL-LVH-16 

SKI  1-SKI 14CC(LM) 

SKI2-SKI2+CC(LL) 

400  CONTINUE 

KIX1-SKI1/16. 

KIX2-SKI2/16. 

R1-D1/24. 

R2-D2/24. 

TH1X-TH1/12. 

TH2X-TH2/12. 

IF(ICALB  .EQ.  0)  MO-11 

C CALCULATIONS  OF  THE  HEAT  LOSSES  FROM  THE  UNDERGROUND  PIPES 
CALL  PIPEHL(T,R1 .R2.TH1X.TH2X.KIX1 .KIX2.MO.OTX) 

HLOSS-QTX 

IF(ITER.EO.I)  HLOSX-0. 

C DETERMINE  IF  PIPE  HEAT  LOSS  VALUE  HAS  CONVERGED.  OR  CONTINUE 
C ITERATIONS  IF  REQUIRED 

D E LQT-ABS ( H LOSS-H  LOSX ) /HLOSS 

IF(DELQT.LE.  0.010)  GO  TO  2010 

ITER-ITER+1 

HLOSX-HLOSS 

GO  TO  38 

2000  WRITE  (7.2005) 

2005  FORMAT  (IX, 'THERE  ARE  SOME  ERRORS  IN  INPUT  DATA’) 

2010  IF(ICALB  .EQ.  1)  THEN 
WRITE  (7,185) 

DO  2020  1-1 ,NE 

WRITE  (7.187)  I.(NODE(I,J).J-1 .3).MAT(I).CC(I) 

2020  CONTINUE 
END  IF 

WRITE  (7.330) 

WRITE  (7,285)  (T( I ) , 1-1 .NN) 

CALL  PIPEHL(T,R1 .R2.TH1X.TH2X.KIX1 .KIX2.7.QTX) 

STOP 

END 

SUBROUTINE  TGO(TGX ,PI ,Y) 

C THIS  SUBROUTINE  CALCULATES  THE  UNDISTURBED  EARTH  TEMPERATURES 
C AT  VARIOUS  DEPTHS 

DIMENSION  TGX(12.6),Y(150) 

READING  IN  THE  ANUAL  AVERAGE  TEMPERATURE  AND  AMPLITUDE  OF  THE 
MONTHLY  NORMAL  TEMPERATURE  CYCLE  OF  THE  SITE.  IN  DEG  F,  AND 
THERMAL  DIFFUSIVITY  OF  SOIL,  IN  FT**2/H. 

READ  (8.*)  AO.BO.DIFF 
W-2. •PI/12. 

WZ-2 . •PI/(8760*DI FF*2) 

ZZ-SQRT(WZ) 

DO  1 1-1,12 
DO  1 J-1 .6 
Z-ZZ*Y(31-J) 

1 TGX(I.J)-AO+BO*EXP(-Z)*SIN(W*(I-3)-Z) 

RETURN 

END 

SUBROUTINE  TGXX(T.TGX, MONTH) 
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C THIS  SUBROUTINE  PROVIDES  OUTER  BOUNDARY  TEMPERATURES  OF  EARTH  REGION 
DIMENSION  T(150) .TGX(12.6) 

T(30)-TGX(MONTH.1) 

DO  1 I«1 .8 
II-I+30 

1 T(II)-=T(30) 

DO  5 1-2,6 
I15-I+15 
JI-31-I 

T(I15)-TGX(MONTH. I) 

T(JI)-TGX(MONTH.I) 

5 CONTINUE 

DO  10  1-1 .3 
I21-I+21 

10  T(I21 )-TGX(MONTH.6) 

RETURN 

END 

SUBROUTINE  INSULK(TM.C) 

C THIS  SUBROUTINE  DETERMINES  THE  THERMAL  CONDUCTIVITY  OF  PIPE 
C INSULATION  (CALCIUM  SILICATE)  AS  A FUNCTION  OF  THE  MEAN 
C TEMPERATURE. 

REAL  KN(16).KINS 
DIMENSION  TN(16) 

DATA  KN  /0. 375, 0.40, 0.42, 0.45, 0.48, 0.50, 0.53, 0.555, 0.58,0. 61 . 
k 0.63,0.66,0.68,0.74,0.82,0.90/ 

DO  5 J-1,16 
IF(J  .LE.  13)  THEN 

TN(J)-100.+(J-1 )»50. 

ELSE 

TN(J)-700.+(J-13)*100. 

END  IF 

5 CONTINUE 

IF(TM  .GT.  TN(1))  GO  TO  10 
KINS-KN(I) 

GO  TO  100 

10  IF(TM  .LT.  TN(16))  GO  TO  20 

KINS=KN(16) 

GO  TO  100 
20  DO  50  1-1,15 

T1-=TM-TN(I) 

IF(T1  .NE.  0.)  GO  TO  30 
KINS=KN(I) 

GO  TO  100 

30  T2-TN(I+1)-TM 

IF(T2  .NE.  0.)  GO  TO  40 
KINS-KN(I+1) 

GO  TO  100 
40  P-T1*T2 

IF(P  .LT.  0.)  GO  TO  50 

KINS-KN(I)+T1*(KN(I+1)-KN(I))/(TN(I+1)-TN(I)) 

GO  TO  100 
50  CONTINUE 

100  C-KINS/12. 

RETURN 

END 

SUBROUTINE  SOILK(TM,KG,C) 

C THIS  ROUTINE  DETERMINES  THE  THERMAL  CONDUCTIVITY  OF  SOIL  AS  A 
C FUNCTION  OF  MEAN  TEMPERATURES. 

REAL  K(14),KG 
DIMENSION  TX(14) 

DATA  K/1 .1 ,1 . 1 ,1 .1 ,1 .0,0.4,0.31 ,0.25,0.19,0.15,0. 11 ,0.09,0.07, 
k 0.05,0.05/ 

DO  1 1-1,14 

1 TX(I)-50.+(I-1)*25. 

IF(TM.GT.TX(1))  GO  TO  5 
ZK-1 . 1 
GO  TO  50 

5 IF(TM.LT.TX(14))  GO  TO  10 

ZK-0.05 
GO  TO  50 

10  DO  25  1-1,13 
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T1-TM-TX(n 

IF(T1 .NE.0)  GO  TO  15 

ZK-K(I) 

GO  TO  50 
15  CONTINUE 

T2-=TM-TX(I+1) 

IF(T2.NE.0.)  GO  TO  20 
ZK-K(I+1) 

GO  TO  50 
20  CONTINUE 

P*T1*T2 

IF(P.GT.0)  GO  TO  25 
ZK-K(I+1)+T2*(K(I+1)-K(I))/25. 

GO  TO  50 
25  CONTINUE 

50  OZK*KG/(1  .1*12.) 

RETURN 

END 

SUBROUTINE  PIPEHL(T.R1 .R2.TH1 ,TH2.ZKS1 .ZKS2.MO.OT) 

C THIS  SUBROUTINE  CALCULATES  THE  AVERAGE  TEMPERATURE  DROPS  ACROSS  THE 
C PIPE  INSULATIONS  AND  THE  RATES  OF  HEAT  LOSS  FROM  THE  UNDERGROUND 
C PIPES  IN  DIRECTLY  BURIED  CONDUIT  SYSTEM 
DIMENSION  T(150) 

PI-=4.  •ATANCI  . ) 

SUM 1-0. 

SUM2-0 . 

N1-8 

DO  1 1-1  ,N1 
K1-I 
K2-I+8 
K3-I+80 
K4-I+88 

SUM1»SUM1+T(K1 )-T(K3) 

SUM2-SUM2+T ( K2 )-T (K4 ) 

1 CONTINUE 

T1-SUM1/N1 
T2-SUM2/N1 
ZKIS1-ZKS1*12. 

ZKIS2-ZKS2*12. 

Q1-ZKS1*2.*PI*T1/L0G((R1+TH1)/R1) 

Q2-ZKS2*2 . *PI •T2/L0G( (R2+TH2)/R2) 

QT-Q1+Q2 

IF(MO  .EG.  11)  GO  TO  50 
WRITE(M0.5)  ZKIS1.ZKIS2 

5 FORMAT(/'  AVERAGE  VALUES  OF  PIPE  INSULATION  THERMAL*. 

4*  CONDUCTIVITY  KI1  - ’.FID. 3,*  KI2  - *.F10.3. 

k * BTU-IN/H-FT**2-DEG  F *) 

WRITE(MO.10)  T1,T2 

10  FORMAT(/’  AVERAGE  TEMPERATURE  DROPS  ACROSS  INSULATION  : * ,/, 
k * T1-  '.F10.2.*  T2-  •,F10.2.*  DEG  F*) 

WRITE(MO.20)  Q1.Q2,QT 

20  FORMAT (/.2X. 'HEAT  LOSSES  FROM  UNDERGROUND  PIPES  :*/*  Q1-’ , 

4 F10.2.  * Q2-*.F10.2.*  QT-*.F10.2.*  BTU/H-FT*) 

50  RETURN 

END 

SUBROUT INE  PIPEO(X.Y.TRTK.D.F.INXK) 

C THIS  SUBROUTINE  READS  IN  THE  INPUT  DATA  TO  BE  USED  FOR  CALCULATIONS 
C OF  THE  HEAT  LOSSES  FROM  THE  UNDERGROUND  PIPES  AND  GENERATES  X AND  Y- 
C COORDINATES  OF  NODAL  POINTS  FOR  THE  TWO  PIPE  SYSTEM. 

REAL  KII .KIG.KI .KG 
DIMENSION  X(150).Y(150) 

COMMON  /PP/JP\ .TP2.KII .KIG.DI1 .DI2.THI1 .THI2.B1 .B2.S1 .S2.TG. 

4 WW.HY. MONTH 

COMMON  /EKC/D 1 P . D2P . D I C . THK 1 . THK2 
C READ  TEMPERATURE  OF  PIPE  NUMBERS  1 AND  2,  IN  DEG  F 
READ  (8.*)  TP1.TP2 

C READ  THERMAL  CONDUCTIVITY  OF  THERMAL  INSULATION  AND  SOIL. 

C RESPECTIVELY.  IN  BTU-IN./H-FT**2  - DEG  F.  AND  INDEX  OF  THERMAL 
C CONDUCTIVITY  : INXK  » 0 CONSTANT  THERMAL  CONDUCTIVITY 
C - 1 TEMPERATURE  DEPENDENT  THERMAL  CONDUCTIVITY 

READ  (8.*)  KII. KI6. INXK 
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C READING  IN  THE  OUTSIDE  DIAMETERS  OF  STEEL  PIPES  1 AND  2.  IN  INCHES 
READ  (8.*)  Dll .DI2 

C READING  IN  THE  THICKNESS  OF  THERMAL  INSULATION  USED  FOR  PIPES  1 
C AND  2.  RESPECTIVELY.  IN  INCHES 
READ  (8.*)  THI1 .THI2 

C READ  IN  THE  INSIDE  DIAMETER  AND  THICKNESS  OF  OUTER  CASING.  AND  THE 
C DEPTH  OF  ITS  CENTER  BELOW  GROUND  SURFACE.  IN  INCHES. 

READ  (8.*)  DIAC.THKC.DEPC 

C READING  IN  THE  VERTICAL  DISTANCES  (IN  FT.)  FROM  HORIZONTAL 

C CENTERLINE  OF  THE  INNER  EARTH  REGION  TO  CENTERS  OF  PIPES  1 AND 

C 2.  RESPECTIVELY.  AND  THE  AVERAGE  EARTH  TEMPERATURE.  IN  DEG  F. 

READ  (8.*)  S1.S2.TG 

C READ  IN  THE  THICKNESS  OF  EARTH  COVER. WIDTH  AND  DEPTH  OF  OUTER 

C EARTH  REGION  SURROUNDING  THE  UNDERGROUND  SYSTEM.  IN  FT. 

READ  (8.*)  E.WW.HY 
WRITE(7.10)  TP1 .TP2.KII .KIG.DI1 .DI2 
10  FORMATC  TP1  TP2  KI  KG  D1  D2’/6F7.2) 
WRITE(7.20)  THI1 .THI2. DIAC.THKC.DEPC. SI .S2.TG 
20  FORMATC  THI1  THI2  DIAC  THKC  DEPC  SI  S2  TG’. 
ft  /8F7.2) 

WRITE(7.30)  E.WW.HY. MONTH 

30  FORMATC  E WW  HY  MONTH  V3F7.2.I7) 

C READ  IN  THE  INSIDE  WIDTH  AND  HEIGHT  OF  THE  INNER  EARTH  REGION. 

C FOR  LOOSE-FILL  INSULATION  SYSTEM.  IN  FT. 

READ  (8.*)  A.B 

C CHANGE  TO  ENGINEERING  UNITS 
D1-DI1/12. 

R1«=D1*0.5 

D2-=DI2/12. 

R2-D2*0.5 

D1P-DI1/12. 

D2P-DI2/12. 

DEP-DEPC/12. 

DIODIAC/12. 

KI*KII/12. 

KG=KIG/12. 

W-=A+2*TRTK/12 

H=B+E+D+F 

WRITE(7.40)  W.H.D.F.A.B.WW.HY 

40  FORMATC  W H D F A B WW  HY* 

ft./.8F7.2) 

PI«4.*ATAN(1  .) 

TH1-:THI1/12. 

TH2=THI2/12. 

THCA=THKC/12. 

THK1-THI1/12. 

THK2-THI2/12. 

C DETERMINE  THE  X AND  Y-COORDINATES  OF  EARTH  COVER,  SIDES  AND  BOTTOM 

C BED  OF  THE  INNER  EARTH  REGION  (NODAL  POINTS  39  TO  70) 

DO  50  I-=1  .5.2 

1 65- 1+65 

X(I65)-W-(I-1)*W/4. 

50  Yn65)-E 

DO  60  1-1. 3. 2 

1 66- 1+66 

X(I66)-(W+A)*0.5  - A*(I-1)*0.5 
60  Y(I66)-E 

DO  65  1-1  .4 
I38-I+38 
I 42- I +42 
146-1+46 
I 50- I +50 
I1-I-1 

X(I38)-(W-A)*0.5 

Y(l38)-(D+E)+I1*B/4. 

X(I42)-0.5*(W-A)+I1*A/4. 

Y(I42)-D+E+B 
X(I46)-(W+A)*0.5 
Yn46)-(D+E+B)-1 1 •B/4. 
xn50)-0.5*(W+A)-I1*A/4. 

Y(I50)-D+E 
65  CONTINUE 

DO  70  1-1 .3 
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1 54- 1 +54 
I 58- I +58 
I 62- I +62 
X(I54)-0.0 

Y( I54)-f D+E)+B* ( 1-1 )/2 . 

X(l58)-(W-A)*0.5+(I-1)*A*0.5 

Y(I58)-H 

X( I62)-W 

Y(l62)-(D+E)+(3-I)*B/2. 

70  CONTINUE 

X(58)-0.0 

Yf58)-H 

X(62)-W 

Y(62)-H 

C THE  X AND  Y-COORDINATES  OF  OUTER  BOUNDARY  EARTH  SURROUNDING  THE 
C DIRECTLY  BURIED  CONDUITS  (NODAL  POINTS  17  TO  38. AND  71  TO  80) 

DO  72  1-1 .2 
I16-I+16 
127=1+27 
X(I16)— WW 

Y(I16)=E+(H-E)*(I-1)*0.5 

X(I27)-W+WW 

Y(I27)=E+(H-E)*(2-I)*0.5 
72  CONTINUE 

DO  75  1-1 .3 
118=1+18 
121=1+21 
I24-I+24 
X(I18)— WW 

y(i18)-H+HY*(I-1 )*0.5 
x(l2l)=W*(I-1)*0.5 
y(I21 )»H+HY 
xn24)=w+ww 
Y(I24)»H+HY*(3-I)*0.5 
75  CONTINUE 

DO  77  1-1 .5.2 

131- I+31 

X(I31)»W-(I-1)*W/4. 

77  YU31)=0.0 
DO  78  1=1 .3.2 

1 32-  I +32 

X( I32)=(W+A) *0 . 5-A* ( 1-1 ) *0 . 5 

78  Y(I32)=0.0 
DO  80  1-1 .2 

I 29-1+29 
I36-I+36 

X( I29)-W+WW* (3-1 ) *0 . 5 
Y(l29)-0. 

Xn36)— WW*0.5*I 
Y(I36)-0. 

80  CONTINUE 

DO  82  1-1 .3 
I 70= I +70 
I 77-1+77 
X(I70)— WW*0.5 
Y(I70)-E+(H-E)*(I-1)*0.5 
X(I77)-^¥+WW*0.5 
Y(I77)-E+(H-E)*(3-I)*0.5 
82  CONTINUE 

DO  85  1-1.4 
I73-I+73 

X(I73)-W*(I-1)/3. 

85  Y(I73)-H+HY*0.5 

C X AND  Y-COORDINATES  OF  THE  CENTERS  OF  THE  PIPES 
XC1i^*0.5 
B1-DEP+S1 
YC1-B1 
XC2-W*0.5 
B2-DEP-S2 
YC2-B2 

WRITE(7.90)  XC1 .YC1 .XC2.YC2.B1 .B2 
90  FORMAT  ('  XC1  YC1  XC2  YC2  B1  B2V6F7.3) 

C THE  X AND  Y-COORDINATES  OF  NODAL  POINTS  AT  THE  INNER  AND  OUTER 
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C SURFACES  OF  PIPE  INSULATION  (NODAL  POINTS  1 TO  16.  AND  81  TO 
C 96) 

DO  95  1-1 ,8 

THETA-2. •PI *1/8. 

1 8- 1+8 
180=1+80 
188=1+88 

Xf n-XC1+0.5*D1*SIN(THETA) 

Y(I)-YC1+0.5*D1*COS(THETA) 

X(I8)-XC2+0.5*D2*SIN(THETA) 

Y (1 8 ) =YC2+0 . 5 • D2 • COS ( THET A ) 


X( 

180' 

)=XC1+( 

'TH1+R1 ' 

)*SINI 

f THETA 

Y| 

180; 

l=YC1+( 

TH1+R1 

)•COS( 

THETA 

X( 

188' 

l=XC2+( 

TH2+R2 

I^SINI 

THETA 

Y( 

188 

l=YC2+( 

;TH2+R2, 

(•COS! 

;theta 

CONTINUE 

C THE  X AND  Y-COORDINATES  OF  NODAL  POINTS  AT  THE  INNER  AND  OUTER 
C SURFACES  OF  PIPE  CASINGS  (NODAL  POINTS  97  TO  112) 

DO  97  1=1 .8 

THETA=2.*PI*I/8. 

196=1+96 

1104=1+104 

X(I96)=XC1+DIC*0.5*SIN(THETA) 

Yn96^DEP+DIC*0.5*COS(THETA) 
xhl04)=XC1  + (THCA+DlC*0.5)*SIN(THETA) 
Y(I104)=DEP+(THCA+DIC*0.5)*COS(THETA) 

97  CONTINUE 

C THE  X AND  Y-COORDINATES  OF  NODAL  POINTS  IN  BACK-FILL  SOIL.  OR 

C POURED- IN  INSULATION  SURROUNDING  THE  PIPES  (NODAL  POINTS  113  TO 

C 124) 

YUP-0.5*(Y(108)-Y(53)) 

XLT-0.5*(X(110)-X(4l)) 

YLO=0.5*(Y(45)-Y(112)) 

XRT-0.5*(X(49)-X(106)) 

DO  100  1-1,3 
1112=1+112 
1115=1+115 
1118=1+118 
1121=1+121 

X ( 1 1 1 2 )-0 . 5* ( W+A)-0 . 25*A* I 
Yni12)=(D+E)+YUP 
xni15)=0.5*(W-A)+XLT 
Y(I115)=(D+E)+0.25*B*I 
X( 1 1 1 8)-0 . 5* (W-A)+0 . 25*A* I 
Yni18)=(D+E+B)-YLO 
X ( 1 1 21 )=0 . 5* (W+A)-XRT 
Y(I121)=(D+E+B)-0.25*B*I 
100  CONTINUE 
RETURN 
END 

SUBROUTINE  TWDPIP( IREPT) 

C THIS  SUBROUTINE  DETERMINES  THE  HEAT  LOSSES  FROM  TWO  PIPES  TO  THE 
C UNDERGROUND  SURROUNDING  THE  HEAT  DISTRIBUTION  SYSTEM. 

REAL  KII.KIG 

COMMON  /PP/T1 .T2.KII .KIG.DI1 .DI2.THI1 ,THI2.B1 ,B2,S1 .S2.TG. 

* WW.HY. MONTH 
PI-4.*ATAN(1 .) 

X1=2.*PI 

R1-DI1/24. 

R2-DI2/24. 

TH1X-THI1/12. 

TH2-THI2/12. 

ZK1-KII/12. 

ZK2-ZK1 

D1-B1 

D2-B2 

ZKS-KIG/12. 

DO  10  1-1 .IREPT 
2 TH1-TH1X+0.1*(I-1) 

S-S1+S2 

A-R1+R2+TH1+TH2+0 . 05 
THI1-TH1*12. 
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IF(A  .LT.  S)  A-S 
C1-X1*2K1/L0G((R1+TH1)/R1  ) 

C2-X 1 • ZK2/LCX5  ( ( R2+TH2  )/R2 ) 

P1 1-=1  .+C1/(X1*ZKS)*L0G((2*D1  )/(R1+TH1)) 

P12=C2/fX1*ZKS)*LOGf (A*A+(D1+D2)**2y(A*A+(D1-D2)**2))*0.5 
P21=C1/(X1*ZKS)*LOG((A*A+(D1+D2)**2)/(A*A+(D1-D2)**2))*0.5 
P22=1 .+C2/(X1*ZKS)*L0G((2*D2)/(R2+TH2)) 

DEL=P12*P21-P11*P22 
ZKP1-=C1*(P12-P22)/DEL 
ZKP2=C2* (P21-P1 1 )/DEL 
TP1-=(P12*T2-P22*Tiy(P12-P22} 

TP2=(P21 •T1-P1 1 •T2)/(P21-P1 1 ) 

Q1=ZKP1.(TP1-TG) 

Q2=ZKP2*(TP2-TG) 

0T=Q1+Q2 

TAVG=(T1+T2)*0.5 

ZK-QT/(TAVG-TG) 

WRITE(7.6)  Dll .DI2.S1 .S2.THI1 .KII .KIG.TI ,T2 
6 FORMAT(’  Dll  DI2  SI  S2  THI1  KII  KIG  TP1  TP2*. 
ft/.7F6.2. 1X.2F6.0) 

WRITE(7,8)  Q1 .Q2.QT,ZK 

8 FORMATC  01  02  OT  KP  •/.3F7.2.2X.F6.3/) 

10  CONTINUE 

RETURN 
END 

SUBROUTINE  EOUIKO(TAS.TDEL.KASP) 

C THIS  ROUTINE  CALCULATES  EOUIVALENT  THERMAL  CONDUCTIVITY  OF  AIR 
C SPACE  SURROUNDING  INSULATED  PIPES  INSIDE  THE  OUTER  CASING. 

REAL  KASP 

COM^N  /EKC/D 1 P . D2P . D I C . THK 1 . THK2 
PI«4.*ATAN(1  .) 

C CALCULATE  THERMAL  CONDUCTIVITY.  IN  BTU-FT/H-FT**2-DEG  F.  AND 
C KINEMATIC  VISCOSITY.  IN  FT**2/S.  OF  AIR 
THKAIR-0.01319  + TAS*2.5E  -5 
VAIR=1.2624E  -4  + TAS*5.4E  -7 

C CALCULATE  THE  EFFECTIVE  DIAMETER  OF  THE  INSULATED  PIPES. 

C IN  FT.  AND  THE  CHARACTERISTIC  LENGTH  OF  AIR  SPACE 
DEFPIP*(D1P42 . •THKI )+(D2P+2 . •THK2) 

CL-DIC-DEFPIP 

C CALCULATE  THE  PRANDTL  NUMBER  OF  AIR  . AND  GRASHOF  NUMBER  AND 
C EOUIVALENT  THERMAL  CONDUCTIVITY.  IN  BTU-IN/H-FT**2-DEG  F.  OF 
r AIR  ^PArr 

PRANTL-0.71849  - TAS  • 1.275E  -4 

GRAS0F-=32.2  * TDEL  • (CL**3 . )/( (VAIR**2 . ) • (TAS+459 . 7) ) 

KASP«=1 2 . ATHKAIRtD . 42*  (PRANTL*GRASOF)  ••0.219 

RETURN 

END 

SUBROUTINE  SOLVLE(A.N.NP. INDX, W.B) 

C GIVEN  AN  NXN  MATRIX  A.  WITH  PHYSICAL  DIMENSION  NP,  THIS  ROUTINE 
C REPLACE  IT  BY  THE  LU  DECOMPOSITION  OF  A ROWWISE  PERMUTATION  OF 
C ITSELF.  INDX  IS  AN  OUTPUT  VECTOR  WHICH  RECORD  THE  ROW  PERMUTATION 

C EFFECTED  BY  THE  PARTIAL  PIVOTING;  W IS  VECTOR  OF  SCALING  FACTORS. 

C 

C THIS  ROUTINE  IS  USED  TO  SOLVE  THE  LINEAR  SET  OF  EOUATIONS  : 

C [A][X]-[B] 

C 

DIMENSION  A(NP.NP) . INDX(N) .VV(N) .B(N) 

C 

C FORM  IMPLICIT  SCALING  VECTOR  W 
C 

DO  12  1-1  .N 
AAMAX  - 0.0 
DO  11  J-1.N 

IF(ABS(A(I .J)) .GT. AAMAX)  AAMAX=ABS(A(I . J)) 

11  CONTINUE 

IF(AAMAX.EQ.0.)  THEN 
WRITE(7.100)  I 

100  FORMAT (IX. ‘ERROR: SINGULAR  MATRIX  - ZERO  ROW  : ROW’. 15) 

RETURN 
END  IF 

VV(I)  « 1 .0/AAMAX 
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12  CONTINUE 


CROUT  METHOD:  LOOP  OVER  COLUMNS 

DO  19  J-=1  .N 
DO  14  1-1 . J-1 
SUM  - A(I , J) 

DO  13  K-1 , 1-1 

SUM  - SUM  - A(I .K)*A(K. J) 

13  CONTINUE 
A(I  , J)  - SUM 

14  CONTINUE 

PIVOT  IMPLEMENTATION 

AAMAX  - 0.0D0 
DO  16  I-J,N 
SUM  - A(I . J) 

DO  15  K-1 .J-1 

SUM  - SUM  - A(I .K)*A(K. J) 

15  CONTINUE 

A(I . J)  - SUM 

DUM  - VV(I)*ABS(SUM) 

IF(DUM.GE.AAMAX)  THEN 
IMAX  = I 
AAMAX  - DUM 
END  IF 

16  CONTINUE 

IF(J.NE. IMAX)  THEN 
DO  17  K-1 ,N 

DUM  - A(IMAX,K) 

AHMAX.K)  - A(J.K) 

A(J.K)  - DUM 

17  CONTINUE 
VV(IMAX)  - VV(J) 

END  IF 

INDXrj)  - IMAX 
IF(A(J .J) .EO.0.0)  THEN 
WRITE(7,110)  J 

0 FORMAT OX. ’ERROR:  SINGULAR  MATRIX  - ZERO  " DIAG  " 

RETURN 
END  IF 

IF(J.NE.N)  THEN 
DUM  - 1 .0/A(J. J) 

DO  18  I-J+1 .N 

A(I.J)  - A(I.J)*DUM 

18  CONTINUE 
END  IF 

19  CONTINUE 

FORWARD  SUBSTITUTION 

II  - 0 
DO  22  1-1 .N 
LL  - INDX(I) 

SUM  - B(LL) 

B(LL)  - B(I) 

IF(II.NE.0)  THEN 
DO  21  J-II .1-1 

SUM  - SUM  - A(I . J)*B(J) 

21  CONTINUE 

ELSE  IF(SUM.NE.0.0)  THEN 
II  - I 
END  IF 
B(I)  - SUM 

22  CONTINUE 

BACKWARD  SUBSTITUTION 

DO  24  I-N.1 .-1 
SUM  - B(I) 

IF(I.LT.N)  THEN 


: ROW’ .15) 
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DO  23  J-I+1 ,N 

SUM  -=  SUM  - A(I  .J)*B(J) 

23  CONTINUE 
END  IF 

B(I)  -=  SUM/A(I.I) 

24  CONTINUE 
RETURN 
END 
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PROGRAM  DIRECT2 

C THIS  IS  A MAIN  PROGRAM  FOR  HEAT  LOSS  ANALYSIS  OF  DIRECTLY  BURIED 
C CONDUIT  UNDERGROUND  HEAT  DISTRIBUTION  SYSTEMS  WITH  INSULATED  PIPES 
C IN  SEPARATE  CASINGS.  BASED  ON  THE  FINITE  ELEMENT  METHOD  USING 
C THREE  - NODE  LINEAR  TRIANGULAR  ELEMENTS. 

C SUBROUTINES  CALLED:  PIPES . TGO . SOI LK . INSULK . TGXX . SOLVLE, PIPEHL . TWOPIP . 
C EQUIKC. 

C INPUT  DATA  FILES:  CDATA1  AND  CDATA2 

C OUTPUT  FILE:  COUTPUT 

C X(n:  THE  X-COORDINATE  OF  NODAL  POINT  I . IN  FT 

C Y(I):  THE  Y-COORDINATE  OF  NODAL  POINT  I, IN  FT 

C (NODE(M, I ) . 1=1 .3) : THREE  NODAL  POINTS  OF  ELEMENT  M 
C M ELEMENT  INDEX 

C NE  TOTAL  NUMBER  OF  ELEMENTS 

C NN  TOTAL  NUMBER  OF  NODAL  POINTS 

C MZ  TOTAL  NUMBER  OF  KNOWN  NODAL  TEMPERATURES 

C C THERMAL  CONDUCTIVITY.  BTU- IN/HR/FT* •2/DEG  F 

C L THICKNESS  OF  THE  ELEMENT.  FT 

C T(I):  THE  TEMPERATURE  OF  NODAL  POINT  I.  IN  DEG  F 
REAL  L.KK.KI .KG.KIX1 ,KIX2.KTCT,KASP.KCAS.KBF 
CHARACTER*4  TITLE05) 

DIMENSION  0(150) .T( 150) .X( 150) ,Yf 150), KK( 150. 150) 

DIMENSION  AS(250) .B2IZ(250) .B3IZ(250) .B2JZ(250) .B2K2(250) . 
ft  B3J2(250) .B3KZ(250) 

DIMENSION  CC(250) .TGX(12.6) .00(150) .NODE(250. 3) ,MAT(250) 
DIMENSION  HIJ(250)  .HJK(250)  .HKI  (250KTIJ(250)  .TJK(250)  , 
ft  TKI (250) .HHl J(250) .HHJK(250) .HHKI (250) . IXCB(250) 

DIMENSION  CK(150. 150)  .00050)  .XT (150)  ,INDX(  150), W(  150) 
COtyMON/PP/TP1  .TP2.KI  .KG.D1  .D2.TH1  .TH2.DP1  .DP2.S1  .S2.TG, 
ft  WW.HY. MONTH 

DIMENSION  TAS(2) .TDEL(2) .KASP(2) .THGP(2) .DIC(2)  .TEPS(2)  . 
ft  TSIC(2) 

COK*ADN  /EKC/D 1 P . D2P . THK 1 . THK2 , E I NS . EGAS 
COMMON  /ST/AO. BO, D IFF 
PI«4.*ATAN(1 . ) 

OPEN  (8,FILE=’CDATA1 ’ ) 

OPEN  (7,FILE=*C0UTPUT’  , STATUS- ' NEW ' . FORIi^’ FORMATTED’ ) 

OPEN  (9,FILE=’CDATA2*) 

C READ  IN  THE  TITLE  OF  THE  PROBLEM  TO  BE  ANALYZED 
READ  (8,2.ERR=2000)  TITLE 

2 FORMAT (15A4) 

WRITE  (7,3)  TITLE 

3 FORMAT  OX,15A4) 

C READ  TOTAL  NUMBER  OF  NODAL  POINTS.  TOTAL  NUMBER  OF  TRIANGULAR 
C ELEMENTS.  TOTAL  NUMBER  OF  KNOWN  NODAL  TEMPERATURES.  AND  THE 
C FIRST  ELEMENT  INDEX  OF  PIPE  INSULATION 
READ  (8.*)  NN.NE.MZ.MINS 

C READ  IN  THE  NUMBER  OF  ITERATIONS  TO  ACCOUNT  FOR  THE  TEMPERATURE 
C EFFECT  ON  INSULATION  AND  SOIL  THERMAL  CONDUCTIVITIES 
READ  (8.*)  MREPT 

C SET  THE  UNIT  NUMBER  OF  THE  PRINTER 
M(V7 

C READ  MONTH  OF  INTEREST  AND  THE  INDEX  FOR  FINITE  ELEMENT  GRID  DATA 
C TO  BE  PRINTED  OUT  : ICALB  - 1 PRINT  OUT  NODAL  COORDINATES 
C - 0 NO  PRINT  OUT 

READ  (8.*)  MONTH. ICALB 

C READ  THE  THERMAL  CONDUCTIVITY  (IN  BTU-IN./H-FT**2  - DEG  F) . 

C THICKNESS  (IN  INCHES)  OF  THE  SIDE,  THE  DEPTH  (IN  FT.)  OF  EARTH 
C COVER  (IN  FT.).  AND  THE  THICKNESS  (IN  FT.)  OF  BOTTOM  BED  OF 
C THE  INNER  EARTH  REGION. 

READ  (8.*)  KTCT.TRTK.D.F 

C READ  IN  THE  ESTIMATED  AVERAGE  TEMPERATURE  OF  AIR  INSIDE  THE 

C ANNULUS  BETWEEN  INSULATED  PIPE  AND  OUTER  CASING,  IN  DEG  F, 

C AND  THE  TEMPERATURE  DIFFERENCE  BETWEEN  THE  INSULATED  PIPE 
C SURFACE  TEMPERATURE  AND  THE  INNER  SURFACE  TEMPERATURE  OF  THE 
C CASING.  IN  DEG  F.  FOR  PIPE  NO.  1 AND  2.  RESPECTIVELY. 

READ  (8,*)  (TAS(I).TDELn).I“1  .2) 

C READ  IN  THE  THERMAL  CONDUCTIVITY  (IN  BTU-IN/H-FT**2-DEG  F)  OF 
C THE  PIPE  CASING  AND  THE  CONDUCTIVITY  OF  POURED-IN  INSULATI(3N 
C MATERIAL  OR  BACK-FILL  SOIL  SURROUNDING  THE  PIPES  IN  THE 
C INNERMOST  REGION. 

READ  (8.*)  KCAS.KBF 

C READ  IN  TOTAL  EWWISSIVITY  OF  OUTER  SURFACE  OF  PIPE  INSULATION 
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C AND  INNER  SURFACE  OF  CONDUIT  CASING.  RESPECTIVELY. AND  THE 
C ESTIMATED  AVERAGE  VALUE  OF  OUTER  SURFACE  TEMPERATURES  FOR 
C PIPES  NO.  1 AND  2.  RESPECTIVELY.  AND  THE  ESTIMATED  INNER 
C SURFACE  TEMPERATURE  OF  THE  CONDUIT  CASING  FOR  PIPES  NO.  1 
C AND  2.  RESPECTIVELY 

READ(8.*)  EINS.ECAS.(TEPS(I).I-1  . 2) . (TSIC(  I ) . I-=1 . 2) 

C READING  IN  INPUT  DATA  FOR  CALCULATIONS  OF  PIPE  HEAT  LOSS  AND 
C GENERATION  OF  THE  COORDINATES  OF  NODAL  POINTS 
CALL  PIPES(X.Y.TRTK.D.F. INXK.THGP.DIC. ITYPE) 

CALL  TWOPIP(I) 

CALL  EOU I KC ( TAS . TDE L . THGP . D I C . T EPS . TS I C . KASP ) 

IF(ICALB.EQ.I)  THEN 
WRITE(7.5) 

5 FORMATC  X(M)  .M=1  .NN') 

WRITE(7.7)  (X(I) . I»1 .NN) 

7 FORMAT (10F7. 2) 

WRITE(7.10) 

10  FORMATC  Y(M)  .M=1  .NN’) 

WRITE(7.7)  (Y(I) .1-1 .NN) 

END  IF 

C CALCULATIONS  OF  UNDISTURBED  EARTH  TEMPERATURES  AT  VARIOUS  DEPTHS 
CALL  TGO(TGX.PI .Y) 

C INITIALIZATION  OF  THE  INDEX  OF  CONVECTION  BOUNDARY  FOR  ELEMENT  N 
DO  12  N-1 .NE 
12  IXCB(N)-0 

C PERFORM  ITERATIONS  TO  ACCOUNT  FOR  THE  TEMPERATURE  EFFECTS  ON  SOIL 
C AND  INSULATION  THERMAL  CONDUCTIVITIES 
DO  24  1-1 .NN 
24  T(I)-TG 

DO  26  1-1 .NE 
HIJI 
HJK( 

HKII 
TIJ( 

TJK( 

TKIi 
HHIJ( 

HHJK( 

HHKK 

26  CONTINUE 

READING  IN  THE  ELEMENT  NUMBER  AND 
MATERIAL  TYPE.  WHICH  INCLUDES 


I)-0. 

-0. 


ITS  NODAL  POINTS  AND  THE 


MAT(J)  - 1 SOIL  IN  INNER  EARTH  REGION 

- 2 PIPE  INSULATION 

- 3 AIR  SPACE  SURROUNDING  THE  PIPE  IN  CASING 

- 4 OUTER  CASING  OF  THE  INSULATED  PIPE 

- 5 POURED-IN  INSULATION  MATERIAL  IN  BACK-FILL 

REGION 

- 6 NATIVE  SOIL  IN  THE  OUTER  EARTH  REGION 


DO  30  1-1. NE 
READ(9.*) 


30 

C 

C 

C 

C 


35 


38 


IF  1 

[MAT! 

;j; 

I.EQ.1) 

IF  ( 

MAT( 

J 

I.EQ.2) 

IF  ( 

MAT( 

j; 

I.EQ.3) 

IF  ({ 

J, 

,GT.122 

IF  (I 
END  IF 

.GT.138 

IF  ( 

;mat( 

».E0.4) 

IF  ( 

MATI 

J 

».EQ.5) 

IF  (MAT! 
CONTINUE 

[j! 

).EQ.6) 

J . (NODE(J .K) .K-1 .3) .MAT(J) 
CC(J)-KTCT/12. 
CC(J)-KI/12. 


.fJ.LT.139))  CC(J)-KASP(1)/12. 
.(J.LT.155))  CCO)-KASP(2)/12. 

-KCAS/12. 

-KBF/12. 

-KG/12. 

JTINUE 

READ  IN  TOTAL  NUMBER  OF  ELEMENTS  HAVING  BOUNDARY  SEGMENTS  SUBJECT 
TO  CONVECTIVE  HEAT  TRANSFER 
READ  (9.*)  NECB 

READ  IN  ELEMENT  NUMBER.  CONVECTIVE  HEAT  TRANSFER  COEFFICIENTS, 

AND  AMBIENT  TEMPERATURES  FOR  THREE  BOUNDARY  SEGMENTS 
DO  35  1-1. NECB 

READ  (9.*)  M.HIJ(M).HJK(M).HKI(M).TIJ(M).TJK(M).TKI(M) 
IXCB(M)-1 
CONTINUE 
ITER-1 

DO  40  1-1, NN 
DO  40  J-1.NN 
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Q(I)-0. 
KKH  . J)-0. 
QQ(I)-0. 
CKfl . J)-0. 

DO(n-0. 

W(I)-1  .0 
INDXU)-! 

40  CONTINUE 


C 

C 
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C 


L-1  . 

DO  180  ,NE 
I-NODEfM, 1 ) 

J-NODE(M.2) 

K-N0DE(M,3) 

IF(MAT(M) .EQ.3)  THEN 

IFffM.GT. 122} .AND. (M. LT. 

IF((M.GT. 138) .AND. (M.LT. 

END  IF 
OCCfM) 

IF  ((INXK.EO.0).OR.(ITER.EQ 
DETERMINE  SOIL  AND  INSULATION  THERMAL  CONDUCTIVITIES  BASED  ON  THE 
MEAN  TEMPERATURES 

TM=(T(I)+T(J)+T(K))/3. 


)«KASPf- 

l)/12 

))  cc(m; 

)«KASP(: 

2)/12 

) GO  TO 

60 

MAT! 

fM) 

.E0.2) 

CALL 

MAT! 

M) 

.EQ.6) 

CALL 

MAT! 

;m) 

.E0.1) 

CALL 

CC(M)-C 


lie 

!.C) 


XI-X 

XJ-=X 

XK-:X 

YI-=Y 

YJ-=Y 

YK=Y 

CXX-C 

CXY-0 . 


CYX-=0. 

CYY-C 

B2I-:YJ-YK 

B3I«XK-XJ 

B2J-YK-YI 

B3J-XI-XK 

B2K=YI-YJ 

B3K*XJ— XI 

CALCULATE  THE  ELEMENT  AREA 
SA-0 . 5* (XJ •B2J+XI •B2 I+XK*B2K) 

SA-ABS(SA) 

A2-SA*2. 

AS(M)-A2 

B2I-B2I/A2 

B3I-B3I/A2 

B2J-B2J/A2 

B3J-B3J/A2 

B2K-B2K/A2 

B3K-B3K/A2 

B2IZ(M)-B2I 

B3IZ(M)-B3I 

B2JZfM)-B2J 

B3JZfM)-B3J 

B2KZ(M)-B2K 

B3KZ(M)-B3K 

BII-SA*L*(B2I*B2I*CXX+B2I*B3I*  CXY+B3 1 • B2 1 *CYX+B3 1 * B3 1 •CYY ) 
BIJ-:SA.L*fB2I*B2J*CXX+B2I*B3J*CXY+B3I*B2J*CYX+B3I*B3J*CYY) 
B IK-SA*  L* ( B2I •B2K*CXX+B2 1 •B3K*CXY+B3I •B2K*CYX+B3I •B3K*CYY) 
B J I -SA*  L»  f B2 J •BZ I •CXX+BZ J •B3 1 •CXY+BSJ •BZ I *CYX+B3 J •B3 1 •CYY) 
B J J«SA*  L* ( B2 J *62 J 'CXX+BZ J •B3 J •CXY+B3 J •BZ J •CYX+BS J •BS J •CYY ) 
BJK-SA^L^(B2J^B2K^CXX+B2J^B3K^CXY+B3J^B2K^CYX+B3J^B3K^CYY) 
BKI«SA^L^(B2K^B2NCXX+B2K^B3NCXY+B3K^B2NCYX+B3K^B3I^CYY) 
BKJ-:SA^L^(B2K^B2J^CXX+B2K^B3J^CXY+B3K^B2J^CYX+B3K^B3J^CYY) 
BKK-SA^L^(B2K^B2K^CXX+B2K^B3K^CXY+B3K^B2K^CYX+B3K^B3K^CYY) 


KK( 

[I.I; 

l-KK( 

.i; 

l+BII 

KK( 

1 .J 

l-KK< 

I .J 

l+BIJ 

KK( 

1 .K 

l-KK( 

I .K 

l+BIK 

KK( 

J.I 

l-KK( 

J.I 

l+BJI 

KK( 

kJ.j' 

l-KK( 

’j.j' 

I+6JJ 
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KKI 

(j.k; 

I-KK( 

;j.K' 

KBJK 

KK| 

K.I 

l-KK( 

K.I 

)+BKI 

KK( 

K.  J 

l=KK( 

K.  J 

l+BKJ 

KK( 

K.K 

l=KK( 

K.K 

l+BKK 

IF( 

;IXCB(M) . 

EQ.0)  GO  TO  130 

C ADDITION  OF  CONVECTION  TERMS  TO  THE  ELEMENT  MATRIX  TO  ACCOUNT 
C FOR  CONVECTION  ON  BOUNDARY 

C READING  IN  CONVECTIVE  HEAT  TRANSFER  COEFFICIENTS  AND  AMBIENT 
C TEMPERATURES  FOR  THREE  BOUNDARY  SEGMENTS 

HHij(M)=Hij(Mi*L*soRTf  (xf  n-xrji).*2+rYrn-Yrj)u*2)/6. 

HHJK(M)«=HJK(M)*L*SQRTf (X( J)-XfK))**2+(Yf j)-Y(K) )**2)/6. 
HHKI(m)=HKI(M)*L*SORT((X(K)-X(I))**2+(Y(K)-YO))**2)/6. 

KK(I . l)=HHIJfM)*2.+HHKI(M)*2.+KK(I .1) 

KKH  , j)«=HHIJfM)+KKf  I . J) 

KKfl ,k)=HHKI fM)+KKf I .K) 

KK( J . l)=HHIJfM)+KK(j. I) 

KKf J . J)=HHIJ(m)*2+HHJK(M)*2.+KK(J. J) 

KKf J .K)-HHJKfM)+KKf J.K) 

KKfK,  I )=hhkhm)+kk(k.  n 
KKfK, j)=HHJKfM)+KK(K.J) 
KK(K.k)=HHJK(M)*2.+HHKI(M)*2.+KK(K.K) 

HHIjiM)=TI J(M)*3. •HHIJfMl 
HHJKfM)<=TJK(M)*3.  •HHJKfM) 

HHKI (M)=TKI (M)*3 . •HHKI (M) 

130  0n)«=0(  n+HHIJ(M)+HHKHM) 

0( J)-q( JI+HHIJfMI+HHJKfM) 

Q(K)=Q(k)+HHJK(M)+HHKI (M) 

180  CONTINUE 

C DETERMINE  IF  FINITE  ELEMENT  GRID  DATA  ARE  TO  BE  PRINTED  OUT 
IF  (nCALB.EQ.I).AND.(IREPT.EQ.MREPT))  THEN 
WRITE(7.185) 

185  FORMATC  M I J K MAT.  C') 

DO  190  1=1 ,NE 

WRITE(7.187)  I . (NODE(I,J).J=1 . 3) ,MAT( I ) ,CC( I ) 

187  FORMAT(1X.5I6.F10.4) 

190  CONTINUE 
END  IF 

C DETERMINE  OUTER  SURFACE  TEMPERATURES  OF  UNDERGROUND  PIPES 
DO  200  1=1.8 
T(I)=TP1 
11=1+8 
T(II)=TP2 
200  CONTINUE 

C DETERMINE  OUTER  BOUNDARY  TEMPERATURES  OF  EARTH  REGION 
CALL  TGXX(T,TGX. MONTH) 

MZ1-MZ+1 
DO  260  I-MZ1.NN 
SUM-0. 

DO  250  J=1 ,MZ 

250  SUM=SUM+KK(I ,J)*T(J) 

OQ(I)-Q(I)-SUM 
260  CONTINUE 

IF(ICALB.EQ.I)  THEN 
WRITE(7.280) 

280  FORMAT (6X. *00  ARRAY’) 

WRITE(7.285)  (QQ(I) , 1=1 .NN) 

285  FORMAT  (5E14.5) 

END  IF 

C RENAMING  OF  MATRICES 
MN-NN-MZ 
DO  300  1-1 ,MN 
K-MZ+I 

DO  290  J=1  ,MN 
KL-MZ+J 

290  CK(I . J)-KK(K.KL) 

XTn)=T(K) 

DQ(I)-QQ(K) 

300  CONTINUE 

C SOLUTION  OF  SIMULTANEOUS  EQUATIONS 
C SET  PHYSICAL  DIMENSION  OF  MATRIX  A 
NP-150 

CALL  SOLVLE(CK.MN.NP.INDX.W.DQ) 

DO  310  1=1 .MN 
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K-M2+I 

T(K)-DO(I) 

310  CONTINUE 

C CALCULATIONS  OF  EQUIVALENT  THERMAL  CONDUCTIVITIES  FOR  AIR 
C SPACES  IN  INSULATED  PIPES 
C 

C CALCULATE  AVERAGE  SURFACE  TEMPERATURE  OF  INSULATED  PIPES 
SU1-0.0 
SU2-0 . 0 
DO  312  1-1 .8 
L1-I+80 
L2-1+88 
SU1-SU1+T(L1 ) 

SU2-SU2+T(L2) 

312  CONTINUE 

TEPS(1)-SU1/8. 

TEPS(2)-SU2/8. 

C CALCULATE  THE  INNER  SURFACE  TEMPERATURE  OF  THE  CASING 
SU3=0 . 0 
SUA-0.0 
DO  314  1-1 .8 
L3-I+96 
L4=I+104 
SU3=SU3+T(L3) 

314  SU4=SU4+T(L4) 

TSICn)=SU3/8. 

TSIC(2)=SU4/8. 

C DETERMINE  THE  TEMPERATURE  DIFFERENCE  BETWEEN  THE  INSULATED  PIPE 
C AND  THE  INNER  SURFACE  OF  THE  CASING  FOR  PIPES  1 AND  2.  RESPECTIVELY. 
DO  316  1-1 .2 

TAS(I)=(TEPS(I)+TSIC(I))/2. 

TDEL( I )-ABS(TEPS( I )-TSIC( I ) ) 

316  CONTINUE 

CALL  EOUIKC(TAS.TDEL.THGP.DIC.TEPS.TSIC.KASP) 

WRITE(7.320)  (I.KASP(I).I-I .2) 

320  FORMATC  KASPf* . 1 1 . * . F10 . 4.2X . ’ (BTU-IN./H-FT**2-DEG  F) ’/. 

ft  * KASP^ . II , * .F10.4.2X. ’ (BTU-IN./H-FT**2-DEG  F) ■ ) 

330  FORMATC  TEMPERATURE  ARRAY  : T(I).  1-1, NN  ') 

C CALCULATE  THE  MEAN  VALUES  OF  INSULATION  THERMAL  CONDUCTIVITY  FOR 
C PIPES  1 AND  2 
350  SKI1-0. 

SKI  2=0. 

DO  400  LN=1 ,16 
LM-MINS+LN-1 
LL-LM+16 

SKI  1-SKI 1+CC(LM) 

SKI2-SKI2+CC(LL) 

400  CONTINUE 

KIX1-SKI1/16. 

KIX2-SKI2/16. 

R1-D1/24. 

R2-D2/24. 

TH1X-TH1/12. 

TH2X-TH2/12. 

IFdCALB  .EQ.  0)  THEN 
MO=11 

IF(IREPT  .EQ.  MREPT)  M07 
END  IF 

C CALCULATIONS  OF  THE  HEAT  LOSSES  FROM  THE  UNDERGROUND  PIPES 
CALL  PIPEHL(T,R1 ,R2,TH1X,TH2X,KIX1 ,KIX2,MO,QTX) 

HLOSS-OTX 

IF(ITER.EO.I)  HLOSX-0. 

C DETERMINE  IF  PIPE  HEAT  LOSS  VALUE  HAS  CONVERGED,  OR  CONTINUE 
C ITERATIONS  IF  REQUIRED 

DE  LQT-ABS ( H LOSS-H  LOSX ) /H  LOSS 

IF(DELQT.LE.  0.010)  GO  TO  2010 

ITER-ITER+1 

HLOSX-HLOSS 

GO  TO  38 

2000  WRITE  (7,2005) 

2005  FORMAT  (IX, ’THERE  ARE  SOME  ERRORS  IN  INPUT  DATA’) 

2010  IF(ICALB  .EQ.  1)  THEN 
WRITE  (7,330) 
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WRITE  (7.285)  (T( I ) . 1-1 . NN) 

END  IF 

STOP 

END 

SUBROUTINE  TGO(TGX.PI .Y) 

C THIS  SUBROUTINE  CALCULATES  THE  UNDISTURBED  EARTH  TEMPERATURES 
C AT  VARIOUS  DEPTHS 

DIMENSION  TGX(12.6) .Y(150) 

READING  IN  THE  ANUAL  AVERAGE  TEMPERATURE  AND  AMPLITUDE  OF  THE 
MONTHLY  NORMAL  TEMPERATURE  CYCLE  OF  THE  SITE.  IN  DEG  F.  AND 
THERMAL  DIFFUSIVITY  OF  SOIL.  IN  FT**2/H. 

READ  (8.*)  AO.BO.DIFF 
W=2. •PI/12. 

WZ=2.*PI/(8760*DIFF*2) 

ZZ=SQRT(WZ) 

DO  1 1-1 .12 
DO  1 J=1 .6 
Z=ZZ*Y(31-J) 

1 TGX(I . J)-AO+BO*EXP(-Z)*SIN(W*(I-3)-Z) 

RETURN 

END 

SUBROUTINE  TGXX(T.TGX. MONTH) 

C THIS  SUBROUTINE  PROVIDES  OUTER  BOUNDARY  TEMPERATURES  OF  EARTH  REGION 
DIMENSION  T(150) .TGX(12.6) 

T(30)=TGX(MONTH. 1) 

DO  1 1=1.8 
II-I+30 

1 T(II)=T(30) 

DO  5 1=2.6 
115=1+15 
JI=31-I 

Tf I15)=TGX(MONTH. I) 

T(JI)-TGX(MONTH.I) 

5 CONTINUE 

DO  10  1=1 .3 
121=1+21 

10  T(I21)=TGX(MONTH.6) 

RETURN 

END 

SUBROUTINE  INSULK(TM.C. ITYPE) 

C THIS  SUBROUTINE  DETERMINES  THE  THERMAL  CONDUCTIVITY  OF  PIPE 
C INSULATIONS  (CALCIUM  SILICATE  AND  MINERAL  WOOL)  AS  A FUNCTION 
C OF  THE  MEAN  TEMPERATURE. 

REAL  KN(16).KINS 
DIMENSION  TN(16) 

DATA  KN  /0. 375. 0.40. 0.42. 0.45. 0.48. 0.50. 0.53. 0.555. 0.58.0. 61 . 
k 0.63.0.66.0.68.0.74.0.82.0.90/ 

C DETERMINE  THE  THERMAL  CONDUCTIVITY  OF  CALCIUM  SILICATE 
IF  (ITYPE  .EQ.  1)  THEN 
DO  5 J=1 .16 
IF(J  .LE.  13)  THEN 
TN(J)=100.+(J-1)*50. 

ELSE 

TN(J)-700.+(J-13)*100. 

END  IF 

5 CONTINUE 

IF(TM  .GT.  TN(1))  GO  TO  10 
KINS-KN(I) 

GO  TO  100 

10  IF(TM  .LT.  TN(16))  GO  TO  20 

KINS-KN(16) 

GO  TO  100 

20  DO  50  1-1.15 

T1«TM-TN(I) 

IF(T1  .NE.  0.)  GO  TO  30 
KINS-KN(I) 

GO  TO  100 

30  T2-TN(I+1)-TM 

IF(T2  .NE.  0.)  GO  TO  40 
KINS-KN(I+1) 
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GO  TO  100 
P»T1*T2 

IF(P  .LT.  0.)  GO  TO  50 

KINS-KN(I)+T1.(KN(I+1)-KN(I))/(TN(I+1)-TN(I)) 

GO  TO  100 
50  CONTINUE 
100  0:KINS/12. 

C DETERMINE  THE  THERMAL  CONDUCTIVITY  OF  MINERAL  WOOL 
ELSE  IF  (ITYPE  .EQ.  2)  THEN 

KINS-0.2420  + TM*1.5501E  -4  4 TM*TM*7 . 5000E  -7 
OK  I NS/1 2. 

ELSE 
RETURN 
END  IF 
RETURN 
END 

SUBROUTINE  SOI LK(TM.KG ,C) 

C THIS  ROUTINE  DETERMINES  THE  THERMAL  CONDUCTIVITY  OF  SOIL  AS  A 
C FUNCTION  OF  MEAN  TEMPERATURES. 

REAL  K(14) ,KG 
DIMENSION  TX(14) 

DATA  K/1 . 1 ,1 .1 .1 . 1 . 1 .0.0.4.0.31 .0.25.0.19,0. 15,0. 11 .0.09,0.07. 
ft  0.05,0.05/ 

DO  1 1=1 .14 

1 TX(I)=50.4(I-1)*25. 

IF(TM.GT.TX(1))  GO  TO  5 
ZK=1 . 1 
GO  TO  50 

5 IF(TM. LT.TX(14))  GO  TO  10 

ZK=0 . 05 
GO  TO  50 

10  DO  25  1=1 .13 
T1=TM-TX(I) 

IF(T1 .NE.0)  GO  TO  15 
ZK=K(I) 

GO  TO  50 
15  CONTINUE 

T2=TM-TX(I41) 

IF(T2.NE.0.)  GO  TO  20 
ZK=K(I41) 

GO  TO  50 
20  CONTINUE 

P=T1*T2 

IF(P.GT.0)  GO  TO  25 

ZK=K( 141 )4T2* (K( I41 )-K( I ) )/25 . 

GO  TO  50 
25  CONTINUE 

50  C=ZK*KG/(1 .1*12. ) 

RETURN 

END 

SUBROUTINE  PIPEHL(T.R1 .R2.TH1 ,TH2.ZKS1 .ZKS2.MO.OT) 

C THIS  SUBROUTINE  CALCULATES  THE  AVERAGE  TEMPERATURE  DROPS  ACROSS  THE 
C PIPE  INSULATIONS  AND  THE  RATES  OF  HEAT  LOSS  FROM  THE  UNDERGROUND 
C PIPES  IN  DIRECTLY  BURIED  CONDUIT  SYSTEM 
DIMENSION  T(150) 

PI-4.*ATAN(1 . ) 

SUM1-0. 

SUM2-0 . 

N1=8 

DO  1 1-1 .N1 
K1  = I 
K2-I48 
K3-I480 
K4-I488 

SUM1-SUMl4TrK1 )-T(K3) 

SUM2-SUM24T ( K2 )-T ( K4 ) 

1 CONTINUE 
T1-SUM1/N1 
T2-SUM2/N1 
ZKIS1=ZKS1*12. 

ZKIS2-ZKS2*12. 
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Q1-ZKS1*2.*PI*T1/L0G((R1+TH1)/R1) 

Q2-ZKS2 • 2 . •P I • T2/L0G ( ( R2+TH2 ) /R2 ) 

0T-Q1+Q2 

IF(MO  .EQ.  11)  GO  TO  50 
WRITE(M0.5)  ZKIS1.ZKIS2 

5 FORMAT(/'  AVERAGE  VALUES  OF  PIPE  INSULATION  THERMAL’. 

ft’  CONDUCTIVITY  KI1  -=  ’.F10.3.’  KI2  -=  ’.F10.3. 

ft  ’ BTU-IN/H-FT**2-DEG  F ’) 

WRITE(MO. 10)  T1 .T2 

10  FORMAT(/’  AVERAGE  TEMPERATURE  DROPS  ACROSS  INSULATION  : ’ ./. 
ft  ’ T1-=  ’.F10.2.’  T2=  ’.F10.2.’  DEG  F’) 

WRITE(MO.20)  01 .Q2.QT 

20  FORMAT (/.2X. ’HEAT  LOSSES  FROM  UNDERGROUND  PIPES  : ’/’  Q1«’ . 

ft  F10.2.  ’ Q2=’.F10.2.’  OT=’.F10.2.’  BTU/H-FT’) 

50  RETURN 

END 

SUBROUTINE  PIPES(X.Y.TRTK.D.F. INXK.THGP.DIC, ITYPE) 

C THIS  SUBROUTINE  READS  IN  THE  INPUT  DATA  TO  BE  USED  FOR  CALCULATIONS 
C OF  THE  HEAT  LOSSES  FROM  THE  UNDERGROUND  PIPES  AND  GENERATES  X AND  Y- 
C COORDINATES  OF  NODAL  POINTS  FOR  THE  TWO  PIPE  SYSTEM. 

REAL  KI I .KIG.KI ,KG 

DIMENSION  X(150) .Y(150) .THGP(2) .DIC(2) 

COlwMON  /PP/TP1  .TP2.KII  .KIG.DII  .DI2.THI1  ,THI2.B1  ,B2,S1  .S2.TG. 
ft  WW.HY, MONTH 

COlwWON  /EKC/D 1 P , D2P , THK 1 . THK2 , E I NS , ECAS 
C READ  TEMPERATURE  OF  PIPE  NUMBERS  1 AND  2.  IN  DEG  F 
READ  (8.*)  TP1 .TP2 

C READ  THERMAL  CONDUCTIVITY  OF  THERMAL  INSULATION  AND  SOIL. 

C RESPECTIVELY.  IN  BTU-IN ./H-FT**2  - DEG  F.  AND  INDEX  OF  THERMAL 
C CONDUCTIVITY  : INXK  = 0 CONSTANT  THERMAL  CONDUCTIVITY 
C = 1 TEMPERATURE  DEPENDENT  THERMAL  CONDUCTIVITY 

C AND  THE  MATERIAL  TYPE  OF  PIPE  INSULATION.  WHICH  INCLUDES 
C ITYPE  = 1 CALCIUM  SILICATE 

C =2  MINERAL  WOOL 

READ  (8.*)  KII .KIG. INXK. ITYPE 

C READING  IN  THE  OUTSIDE  DIAMETERS  OF  STEEL  PIPES  1 AND  2.  IN  INCHES 
READ  (8.*)  DI1.DI2 

C READING  IN  THE  THICKNESS  OF  THERMAL  INSULATION  USED  FOR  PIPES  1 
C AND  2.  RESPECTIVELY.  IN  INCHES 
READ  (8.*)  THI1.THI2 

C READ  IN  THE  THICKNESSES  OF  AIR  GAP  AND  OUTER  CASING  USED  FOR 
C PIPES  1 AND  2.  RESPECTIVELY.  IN  INCHES. 

READ  (8.*)  (THGP(I),I=1 .2).THCA1 .THCA2 
C READ  THE  DEPTHS  FROM  GROUND  SURFACE  TO  THE  CENTERS  OF  PIPES  1 AND 
C 2.  RESPECTIVELY.  IN  FT. 

READ  (8.*)  B1.B2 

C READING  IN  THE  HORIZONTAL  DISTANCES  (IN  FT.)  FROM  VERTICAL 
C CENTERLINE  OF  THE  INNER  EARTH  REGION  TO  CENTERS  OF  PIPES  1 AND 
C 2.  RESPECTIVELY.  AND  THE  AVERAGE  EARTH  TEMPERATURE.  IN  DEG  F. 

READ  (8.*)  S1.S2.TG 

C READ  IN  THE  THICKNESS  OF  EARTH  COVER. WIDTH  AND  DEPTH  OF  OUTER 
C EARTH  REGION  SURROUNDING  THE  UNDERGROUND  SYSTEM.  IN  FT. 

READ  (8.*)  E. WW.HY 
WRITE(7.10)  TP1 .TP2.KII.KIG,DI1.DI2 
10  FORMAT(’  TP1  TP2  Kl  KG  D1  D2’/6F7.2) 
WRITE(7.20)  THIl .THI2.B1 .B2,S1 .S2.TG 
20  FORMAT(’  THIl  THI2  DPI  DP2  SI  S2  TG*. 
ft  /7F7.2) 

WRITE(7.25)  (THGP(I).I-=1  .2).THCA1  .THCA2 
25  FORMAT (’  THGP1  THGP2  THCA1  THCA2  ’/4F7.2) 

WRITE(7.30)  E. WW.HY. MONTH 

30  FORMAT(’  E WW  HY  MONTH  ’/3F7.2.I7) 

C READ  IN  THE  INSIDE  WIDTH  AND  HEIGHT  OF  THE  INNER  EARTH  REGION. 

C FOR  LOOSE-FILL  INSULATION  SYSTEM,  IN  FT. 

READ  (8.*)  A.B 

C CHANGE  TO  ENGINEERING  UNITS 
D1-DI1/12. 

R1-D1*0.5 

D2-DI2/12. 

R2-D2-0.5 

D1P-DI1/12. 

D2P-DI2/12. 
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DICfn-fDI1+2.*THI1+2.*THGP(1))/12. 

DIC(2)-(DI2+2.*THI2+2.*THGP(2))/12. 

KI-KII/12. 

KG-KIG/12. 

W=A+2*TRTK/12 

H^B+E+Ot-F 

WRITE(7.40)  W,H.D,F,A,B.WW.HY 

40  FORMAT  (’W  H D F A BWWHY’ 

4./.8F7.2) 

PI*4. ,ATAN(1 . ) 

TH1=THI1/12. 

TH2=THI2/12. 

THG1=THGP(1 )/12. 

THG2=THGP(2)/12. 

THC1=THCA1/12. 

THC2=THCA2/12. 

THK1=TH11/12. 

THK2=THI2/12. 

C DETERMINE  THE  X AND  Y-COORDINATES  OF  EARTH  COVER.  SIDES  AND  BOTTOM 

C BED  OF  THE  INNER  EARTH  REGION  (NODAL  POINTS  39  TO  70) 

DO  50  1=1 .5.2 
165=1+65 

X(I65)=W-(I-1)*W/4. 

50  Y(I65)=E 

DO  60  1=1 .3.2 
166=1+66 

X(I66)=(W+A)*0.5  - A*(I-1)*0.5 
60  Y(I66)=E 

DO  65  1=1 .4 
138=1+38 
142=1+42 
146=1+46 
150=1+50 
11=1-1 

X(I38)=(W-A)*0.5 
Yf I38)=(D+E)+I1*B/4. 

X( 142 )=0 . 5* (W-A)+I 1 •A/4 . 

Yn42)=D+E+B 
Xn46)=(W+A)*0.5 
Yn46)=(D+E+B)-1 1 •B/4. 

Xn50)=0.5*(W+A)-I1*A/4. 

Y(I50)=D+E 
65  CONTINUE 

DO  70  1=1.3 
154=1+54 
158=1+58 
162=1+62 
X(I54)»0.0 

Yn54)=(D+E)+B*(I-1)/2. 
xn58)=(W-A)*0.5+(I-1)*A*0.5 
Yn58)»H 
Xf I62)-W 

Y(l62)=(l>+E)+(3-I)*B/2. 

70  CONTINUE 

Xf58)=0.0 

Y(58)«H 

Xf62)=W 

Y(62)=H 

C THE  X AND  Y-COORDINATES  OF  OUTER  BOUNDARY  EARTH  SURROUNDING  THE 
C DIRECTLY  BURIED  CONDUITS  (NODAL  POINTS  17  TO  38. AND  71  TO  80) 

DO  72  1=1.2 
116=1+16 
127=1+27 
X(I16)=-WW 

Y(I16)-E+(H-E)*(I-1)*0.5 

xn27)=w+ww 

Y(I27)=E+(H-E)*(2-I)*0.5 
72  CONTINUE 

DO  75  1-1 .3 
118=1+18 
I21-I+21 
I24-I+24 
X(I18)— WW 
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Y(I18)-H+HY*(I-1)*0.5 
Xn21  )-W*(I-1)*0.5 

Yn21  Wh+hy 
xn24)-w+ww 
Y(I24)-H+HY*(3-I)*0.5 
75  CONTINUE 

DO  77  1=1 .5.2 
131=1+31 

X(I31)»W-(I-l)*W/4. 

77  Y(I31)=0.0 
DO  78  1=1 .3.2 

132=1+32 

X(I32)=(W+A)*0.5-A*(I-1)*0.5 

78  Y(I32)=0.0 
DO  80  1=1.2 

129=1+29 

136=1+36 

X ( 1 29 )=W+WW* (3-1 ) *0 . 5 
Yn29)=0. 

X(I36)->-WW*0.5«I 

Y(I36H0- 

80  CONTINUE 

DO  82  1=1 .3 
170=1+70 
177=1+77 
X(I70)-^WW*0.5 
Y(I70)=E+(H-E)*(I-1)*0.5 
X(I77)=W+WW*0.5 
Y(I77)=E+(H-E)*(3-I)*0.5 
82  CONTINUE 

DO  85  1=1.4 
173=1+73 

X(I73)=W*(l-1)/3. 

85  Y(I73)=H+HY*0.5 

C X AND  Y-COORDINATES  OF  THE  CENTERS  OF  THE  PIPES 
XC1=W*0.5  - SI 
YC1=B1 

XC2=W*0.5  + S2 
YC2=B2 

WRITE(7.90)  XC1 .YC1 .XC2.YC2 
90  FORMATC  XC1  YC1  XC2  YC274F7.3) 

C THE  X AND  Y-COORDINATES  OF  NODAL  POINTS  AT  THE  INNER  AND  OUTER 
C SURFACES  OF  PIPE  INSULATION  (NODAL  POINTS  1 TO  16.  AND  81  TO 
C 96) 

DO  95  1-1 .8 

THETA=2.*PI*I/8. 

18=1+8 

180=1+80 

188=1+88 

X( I )=XC1+0 . 5*D1 •SIN(THETA) 

Y ( I )=YC 1 +0 . 5*  D 1 -COS ( THETA ) 

X( I8)=XC2+0 . 5*D2*SIN(THETA) 

YfI8)=YC2+0.5*D2*COS(THETA) 


XI 

180] 

l-XC1+( 

'TH1+R1] 

[•SINI 

'THETA 

Y| 

180 

|-YC1+( 

TH1+R1 

1*0061 

THETA 

XI 

188 

l-XC2+( 

TH2+R2 

»*SIN| 

THETA 

Y| 

188 

|-YC2+( 

;TH2+R2; 

[•COS! 

[THETA 

CONTINUE 

C THE  X AND  Y-COORDINATES  OF  NODAL  POINTS  AT  THE  INNER  AND  OUTER 
C SURFACES  OF  PIPE  CASINGS  (NODAL  POINTS  97  TO  128) 

DO  97  1-1 .8 

THETA-2. -PI *1/8. 

196=1+96 

1104=1+104 

1112=1+112 

1120=1+120 

X(I96)-XC1+(THG1+TH1+R1)*SIN(THETA) 

Y(I96)-YC1+(THG1+TH1+R1)*COS(THETA) 

X (1 1 04 )-XC2+ ( THG2+TH2+R2 ) • S I N ( THETA ) 

Y (1 1 04 )-YC2+ ( THG2+TH2+R2 ) ^COS ( THETA) 

Xfl 1 12)=XC1+(THC1+THG1+TH1+R1 VSINfTHETA) 
Y(I112)-YC1+(THC1+THG1+TH1+R1)*C0S(THETA) 

X U 1 20 ) -XC2+ ( THC2+THG2+TH2+R2 ) • S I N ( TH ET A ) 


73 


Y ( 1 1 20 )-YC2+ ( THC2+THG2+TH2+R2 ) -COS ( THETA ) 

97  CONTINUE 

C THE  X AND  Y-COORDINATES  OF  NODAL  POINTS  IN  BACK-FILL  SOIL.  OR 

C POURED-IN  INSULATION  SURROUNDING  THE  PIPES  (NODAL  POINTS  129  TO 

C 142) 

YUP-0.5*(Y(1 16)-Y(53)) 

YUPR=0.5*(Y(124)-Y(53)) 

IF  (YUPR.  LT.YUP)  YUP*=YUPR 
XLT-=0.5*fX(118)-X(41)) 

YLO=0.5*(Y(45)-Y(120)) 

YLOW=i0.5*(Y(45)-Y(128)) 

IF  (YLOW. LT.YLO)  YLO=YLOW 
XRT=0.5*(X(49)-X(122)) 

DO  100  1=1 .3 
1128=1+128 
1131=1+131 
1134=1+134 
1137=1+137 

X(I128)=0.5*(W+A)-0.25*A*I 
Y(I128)=(D+E)+YUP 
X(I131 )=0.5*(W-A)+XLT 
Y(I131 )=(D+E)+0.25*B*I 
Xni34)=0.5*(W-A)+0.25*A*I 
Y(I134)=(D+E+B)-YLO. 

X(I137)=0.5*(W+A)-XRT 
Y(I137)=(D+E+B)-0.25*B*I 
100  CONTINUE 

DO  120  1=1 ,2 
1140=1+140 
X(I140)=X(53) 

Y(I140)=(D+E+YUP)+(B-YUP-YLO)*I/3. 

120  CONTINUE 
RETURN 
END 

SUBROUTINE  TWOPIP( IREPT) 

C THIS  SUBROUTINE  DETERMINES  THE  HEAT  LOSSES  FROM  TWO  PIPES  TO  THE 
C UNDERGROUND  SURROUNDING  THE  HEAT  DISTRIBUTION  SYSTEM. 

REAL  KII.KIG 

COMMON  /PP/T1 .T2.KII .KIG.DI1 .DI2.THI1 .THI2.B1 .B2.S1 .S2.TG, 

* WW.HY, MONTH 
PIc=4.*ATAN(1  . ) 

X1=2.*PI 

R1=DI1/24. 

R2=DI2/24. 

TH1X=THI1/12. 

TH2=THI2/12. 

ZK1-KII/12. 

ZK2=ZK1 

D1«B1 

D2-B2 

ZKS-KIG/12. 

DO  10  1=1 .IREPT 
2 TH1-TH1X+0. 1*(I-1) 

TH2=TH1 

S=S1+S2 

A-R 1 +R2+TH 1 +TH2+0 . 05 
THI1=TH1*12. 

IF(A  .LT.  S)  A=S 

C1=X1 •ZK1/LOG( (R1+TH1 )/R1 ) 

C2-X1 •ZK2/LOG( (R2+TH2)/R2) 

P1 1 = 1 .+C1/(X1*ZKS)*LOG((2*D1)/(R1+TH1 )) 


P12-C2/( 

;xi*zKs; 

l•LOG( 

(A*A+( 

;di+d2) 

•••2)/ 

^(A*A+( 

;di-d2; 

• ••z) 

)*0.5 

P21=C1/( 

;xi*zKs; 

)•LOG( 

(A*A+( 

:D1+D2j 

•••2)/ 

'(A*A+( 

:di-D2; 

(••Z) 

)*0.5 

P22=1 .+C2/(X1*ZKS)*LOG((2*D2)/(R2+TH2)) 
DEL-P12*P21-P11*P22 
ZKP 1 =C 1 • f P 1 2-P22 ) /DE L 
ZKP2=C2*(P21-P11)/DEL 
TP1=(P12*T2-P22»T1 )/(P12-P22) 

TP2-(P21 •TI-PI 1 •T2)/(P21-P1 1 ) 
01=ZKP1*(TP1-TG) 

02=ZKP2*(TP2-TG) 

OT-Q1+02 
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TAVG-(T1+T2)*0.5 

2K-QT/(TAVG-TG) 

WRITE(7.6)  DI1 .DI2.S1 .S2 .THI 1 .KI I .KIG , T1 .T2 
6 FORMATC  Dll  DI2  SI  S2  THI1  KII  KIG  TP1  TP2’ 
ft/.7F6.2.1X,2F6.0) 

WRITE(7.8)  01 .Q2.0T,ZK 

8 FORMATC  01  02  OT  KP •/.3F7 . 2 .2X , F6 .3/) 

10  CONTINUE 
RETURN 
END 

SUBROUT  I NE  EOU I KC ( TAS . TDEL , THGP , D I C . T EPS , TS I C , KASP ) 

C THIS  ROUTINE  CALCULATES  EOUIVALENT  THERMAL  CONDUCTIVITY  OF  AIR 
C SPACE  IN  THE  ANNULUS  BETWEEN  INSULATED  PIPE  AND  THE  CASING  FOR 
C HORIZONTAL  CONCENTRIC  PIPES. 

REAL  KASP 

DIMENSION  TAS(2) .TDEL(2) ,KASP(2) .THGP(2) ,DIC(2) ,TEPS(2) . 
ft  TSIC(2) ,C0NKA(2) .RADKA(2) ,D0P(2) 

COMMON  /EKC/D 1 P . D2P . THK 1 . THK2 . E I NS . ECAS 
PI=4.*ATANC  .) 

D0P(1}=D1P+2.*THK1 

D0P(2)-=D2P+2.*THK2 

C CALCULATE  THERMAL  CONDUCTIVITY.  IN  BTU-FT/H-FT**2-DEG  F.  AND 
C KINEMATIC  VISCOSITY.  IN  FT*»2/S.  OF  AIR 
DO  30  1=1 .2 

THKAIR=0. 01319  + TAS(I)*2.5E  -5 
VAIR=1.2624E  -4  + TAS(I)*5.4E  -7 
CL=THGP(I)/12. 

C CALCULATE  THE  PRANDTL  NUMBER  OF  AIR  . AND  GRASHOF  NUMBER  AND 
C EOUIVALENT  THERMAL  CONDUCTIVITY.  IN  BTU-IN/H-FT**2-DEG  F,  OF 

C ATR  QPAOP 

PRANTL=0.71849  - TAS(I)*  1 .275E  -4 

GRASOF=32.2  • TDEL( I ) • (CL**3. )/( (VAIR**2 . ) • (TAS( I )+459 . 7) ) 
CONKA( I )=1 2 . •THKAIR*0 . 1 59* (PRANTL*GRASOF) **0 . 272 
C CALCULATE  THE  EMMISIVITY  TERM  AND  EOUIVALENT  THERMAL 
C CONDUCTIVITY  OF  AIRSPACE  DUE  TO  RADIATIVE  TRANSFER.  IN 
C BTU-IN/H-FT**2-DEG  F 

ET0T=1./EINS  + (D0P(I)/DIC{I))*(1 ./ECAS  - 1.) 
TREP=(TEPS(n+459.7)/100. 

TRSI=(TSIC(I)+459.7)/100. 

HRAD=0 .1714/100.* (TREP+TRSI ) * (TREP*TREP+TRS I *TRSI )/ETOT 
RADKA(I)=12.  * HRAD  • CL 
30  CONTINUE 

C CALCULATE  EOUIVALENT  THERMAL  CONDUCTIVITY  OF  AIRSPACE.  IN 

C BTU-IN/H-FT**2-DEG  F 

DO  50  1=1.2 

KASP(I)=CONKA(I)  + RADKA(I) 

WRITE  (7.45)  I.CONKA(I).RADKA(I) 

45  FORMAT!  • I=M5.2X,  * CONKA=’ ,F10.4.2X,  ’ RADKA-’ ,F10.4. 

ft  2X. ’(BTU-IN/H-FT**2-DEG  F)') 

50  CONTINUE 

RETURN 
END 

SUBROUT I NE  SOLVLE ( A . N . NP , I NDX . W . B ) 

C GIVEN  AN  NXN  MATRIX  A.  WITH  PHYSICAL  DIMENSION  NP.  THIS  ROUTINE 
C REPLACE  IT  BY  THE  LU  DECOMPOSITION  OF  A ROWWISE  PERMUTATION  OF 
C ITSELF.  INDX  IS  AN  OUTPUT  VECTOR  WHICH  RECORD  THE  ROW  PERMUTATION 

C EFFECTED  BY  THE  PARTIAL  PIVOTING;  W IS  VECTOR  OF  SCALING  FACTORS 

C 

C THIS  ROUTINE  IS  USED  TO  SOLVE  THE  LINEAR  SET  OF  EQUATIONS  : 

C [A][X]=[B] 

c 

DIMENSION  A(NP.NP) . INDX(N) .VV(N) .B(N) 

C 

C FORM  IMPLICIT  SCALING  VECTOR  W 
C 

DO  12  1=1. N 
AAMAX  = 0.0 
DO  11  J=1 .N 

IF(ABS(A(I ,J)).GT.AAMAX)  AAMAX-ABS(A(I , J)) 

11  CONTINUE 

IF(AAMAX.EQ.0.)  THEN 


75 


ooo  -*  ooo  ooo 


WRITE(7.100)  I 

100  FORMAT( IX, 'ERROR: SINGULAR  MATRIX  - ZERO  ROW  : ROW'. 15) 
RETURN 
END  IF 

W(I)  - 1 .0/AAMAX 

12  CONTINUE 

CROUT  METHOD:  LOOP  OVER  COLUMNS 

DO  19  J=1 .N 
DO  14  I-=1  . J-1 
SUM  -=  A(I  , J) 

DO  13  K=1 . 1-1 

SUM  -=  SUM  - A(I  .K)*A(K.J) 

13  CONTINUE 
A(I  .J)  = SUM 

14  CONTINUE 

PIVOT  IMPLEMENTATION 

AAMAX  -=  0.0D0 
DO  16  I=J.N 
SUM  «=  A(I  .J) 

DO  15  K=1 ,J-1 

SUM  = SUM  - A(I .K)*A(K. J) 

15  CONTINUE 
A(I.J)  = SUM 

DUM  -=  VV(I)*ABS(SUM) 

IF(DUM.GE. AAMAX)  THEN 
IMAX  «=  I 
AAMAX  >=  DUM 
END  IF 

16  CONTINUE 

IF(J .NE. IMAX)  THEN 
DO  17  K=1 .N 

DUM  = A(IMAX.K) 

A(IMAX.K)  •=  A(J.K) 

A(J.K)  = DUM 

17  CONTINUE 
VV(IMAX)  «=  VV(J) 

END  IF 

INDX(J)  = IMAX 
IF(A(J.J).EO.0.0)  THEN 
WRITE(7.110)  J 

0 FORMAT OX. 'ERROR:  SINGULAR  MATRIX  - ZERO  " DIAG  " : ROW', 15) 

RETURN 
END  IF 

IF(J.NE.N)  THEN 
DUM  « 1 .0/A(J.J) 

DO  18  I«J+1 ,N 

A(I.J)  - A(I.J)*DUM 

18  CONTINUE 
END  IF 

19  CONTINUE 

FORWARD  SUBSTITUTION 

II  -=  0 
DO  22  I-=1  .N 
LL  - INDX(I) 

SUM  « B(LL) 

B(LL)  - B(I) 

IF(II.NE.0)  THEN 
DO  21  J-II.I-1 

SUM  - SUM  - A(I .J)*B(J) 

21  CONTINUE 

ELSE  IF(SUM.NE.0.0)  THEN 
II  - I 
END  IF 
B(I)  - SUM 

22  CONTINUE 
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BACKWARD  SUBSTITUTION 

DO  24  I-N. 1 .-1 
SUM  - B(n 
IF(I.LT.N)  THEN 
DO  23  J-I+1 .N 
SUM  = SUM  - A(I , J)*B(J) 

23  CONTINUE 
END  IF 

B(I)  = SUM/A(I.I) 

24  CONTINUE 
RETURN 
END 
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PROGRAM  DBJACKS 

C THIS  IS  A MAIN  PROGRAM  FOR  HEAT  LOSS  ANALYSIS  OF  DIRECTLY  BURIED 
C CONDUIT  UNDERGROUND  HEAT  DISTRIBUTION  SYSTEMS  INSTALLED  AT  FORT 
C JACKSON.  SC,  BASE  ON  THE  FINITE  ELEMENT  METHOD  USING  THREE-NODE 
C LINEAR  TRIANGULAR  ELEMENTS.  THE  UNDERGROUND  SYSTEMS  CONSIST  OF  TWO 
C INSULATED  PIPES  ENCASED  IN  THE  SAME  CONDUIT  CASING. 


SUBROUTINES  CALLED;  PIPEJ .TGO, SOI LK , INSULK. TGXX.SOLVLE.PIPEHL.TWOPIP, 
EOUIKO. 

SDTA1FJ  AND  SDTA2FJ 


INPUT  DATA  FILES: 
OUTPUT  FILE:  SOUTFJ 


Sli; 


THE  X-COORDINATE  OF  NODAL  POINT  I 
THE  Y-COORDINATE  OF  NODAL  POINT  I, 


IN  FT 
IN  FT 


(NODE(M. I) . 1=1 .3) : THREE  NODAL  POINTS  OF  ELEMENT  M 

M ELEMENT  INDEX 

NE  TOTAL  NUMBER  OF  ELEMENTS 

NN  TOTAL  NUMBER  OF  NODAL  POINTS 

MZ  TOTAL  NUMBER  OF  KNOWN  NODAL  TEMPERATURES 

C THERMAL  CONDUCTIVITY.  BTU-IN/HR/FT**2/DEG  F 

L THICKNESS  OF  THE  ELEMENT.  FT 

T(I):  THE  TEMPERATURE  OF  NODAL  POINT  I.  IN  DEG  F 

REAL  L.KK.KI .KG.KIX1 . KIX2 . KTCT .KASP .KCAS .KBF 

CHARACTER*4  TITLE(15) 

DIMENSION  0(150) .T( 150) .X( 150) .YM 50) .KK(150.150) 

DIMENSION  AS(250)  .B2IZU50)  .B3IZ(250)  .B2JZ(250)  .B2KZ(250)  . 
a;  B3JZ(250).B3KZ(250) 

DIMENSION  CC(250) .TGX(12.7) .00(150) .NODE(250.3) .MAT(250) 
DIMENSION  HIJ(250) .HJK(250) .HKI (250) .TI J(250) .TJK(250) . 
k TKI (250) .HHI J(250) .HHJK(250) .HHKI (250) . IXCB(250) 
DIMENSION  CK(150. 150)  .DO(150)  .XT(150)  . INDX(150)  .W(150) 
C0MM0N/PP/TP1 .TP2.KI .KG.D1 .D2.TH1 .TH2.DP1 .DP2.S1 .S2.TG. 
k WW.HY. MONTH 

COM^N  /EKC/D 1 P . D2P . D I C . THK 1 . THK2 
COMMON  /ST/AO. BO. D IFF 
PI=4..ATAN(1 .) 


NP=150 

OPEN  (B.FILE-’SDTAIFJ*) 

OPEN  (7.FILE='S0UTFJ’ .STATUS=*NEW’ .F0RKA=*F0RMATTED’) 

OPEN  (9.FILE=’SDTA2FJ’) 

C READ  IN  THE  TITLE  OF  THE  PROBLEM  TO  BE  ANALYZED 
READ  (8.2.ERR=2000)  TITLE 

2 FORMAT (15A4) 

WRITE  (7.3)  TITLE 

3 FORMAT ( 1 X.1 5 A4) 

C READ  TOTAL  NUMBER  OF  NODAL  POINTS.  TOTAL  NUMBER  OF  TRIANGULAR 
C ELEMENTS.  TOTAL  NUMBER  OF  KNOWN  NODAL  TEMPERATURES.  THE  FIRST 
C ELEMENT  INDEX  OF  PIPE  INSULATION.  THE  FIRST  NODE  INDEXES  OF 
C PIPE  INSULATION  FOR  PIPES  N0.1  AND  NO. 2,  RESPECTIVELY.  AND  THE 
C FIRST  NODE  INDEX  OF  THE  INNER  SURFACE  OF  THE  CONDUIT  CASING. 

READ  (8.*)  NN.NE.MZ.MINS.NINS1 ,NINS2.NSIC 
C SET  THE  UNIT  NUMBER  OF  THE  PRINTER 
MO-7 

C READ  MONTH  OF  INTEREST  AND  THE  INDEX  FOR  FINITE  ELEMENT  GRID  DATA 
C TO  BE  PRINTED  OUT  : ICALB  - 1 PRINT  OUT  INTERMEDIATE  RESULTS 
C DURING  ITERATIONS 

C = 0 NO  PRINT  OUT 


READ  (8.*)  MONTH. ICALB 

C READ  THE  THERMAL  CONDUCTIVITY  (IN  BTU-IN./H-FT**2  - DEG  F) . 

C THICKNESS  (IN  INCHES)  OF  THE  SIDE.  THE  DEPTH  (IN  FT.)  OF  EARTH 
C COVER  (IN  FT.).  AND  THE  THICKNESS  (IN  FT.)  OF  BOTTOM  BED  OF 
C THE  INNER  EARTH  REGION. 

READ  (8.*)  KTCT.TRTK.D.F 

C READ  IN  THE  ESTIMATED  AVERAGE  TEMPERATURE  OF  AIR  INSIDE  THE 
C AIRSPACE  BETWEEN  INSULATED  PIPES  AND  OUTER  CASING.  IN  DEG  F. 

C AND  THE  TEMPERATURE  DIFFERENCE  BETWEEN  THE  INSULATED  PIPE 
C SURFACE  TEMPERATURE  AND  THE  INNER  SURFACE  TEMPERATURE  OF  THE 
C CASING.  IN  DEG  F. 

READ  (8.*)  TAS.TDEL 

C READ  IN  THE  THERMAL  CONDUCTIVITY  (IN  BTU-IN/H-FT**2-DEG  F)  OF 
C THE  PIPE  CASING  AND  THE  CONDUCTIVITY  OF  POURED-IN  INSULATION 
C MATERIAL  OR  BACK-FILL  SOIL  SURROUNDING  THE  PIPES  IN  THE 
C INNERMOST  REGION. 

READ  (8.*)  KCAS. KBF 

C READ  ’N  TOTAL  EMISSIVITIES  OF  OUTER  SURFACE  OF  PIPE  INSULATION 
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AND  INNER  SURFACE  OF  CONDUIT  CASING.  RESPECTIVELY.  AND  THE  ESTIMATED 
AVERGE  VALUES  OF  OUTER  SURFACE  TEMPERATURES  OF  PIPE  INSULATION  AND 
THE  INNER  SURFACE  TEMPERATURES  OF  THE  OUTER  CASING.  IN  DEG  F. 

READ  (8.*)  EINS.ECAS.TEPS.TSIC 

READING  IN  INPUT  DATA  FOR  CALCULATIONS  OF  PIPE  HEAT  LOSS  AND 
GENERATION  OF  THE  COORDINATES  OF  NODAL  POINTS 
CALL  PIPEJ(X.Y.TRTK.D.F. INXK. ITYPE) 

CALL  TWOPIP(I) 

CALL  EOUIKO(TAS.TDEL.EINS. ECAS . T EPS . TS I C . KASP ) 

IF(ICALB.EO.I)  THEN 
WRITE(7.5) 

FORMATC  X(M)  .14=1  .NN') 

WRITE(7.7)  (X(I) . 1=1 .NN) 

FORMAT (10F7. 2) 

WRITE(7. 10) 

FORMATC  Y(M)  .M=1  .NN’ ) 

WRITE(7.7)  (Y(I) . 1=1 .NN) 

END  IF 

CALCULATIONS  OF  UNDISTURBED  EARTH  TEMPERATURES  AT  VARIOUS  DEPTHS 
CALL  TGO(TGX.PI .Y) 

INITIALIZATION  OF  THE  INDEX  OF  CONVECTION  BOUNDARY  FOR  ELEMENT  N 
DO  12  N=1 .NE 
IXCB(N)=0 

PERFORM  ITERATIONS  TO  ACCOUNT  FOR  THE  TEMPERATURE  EFFECTS  ON  SOIL 
AND  INSULATION  THERMAL  CONDUCTIVITIES 
DO  24  1=1 .NN 
T(I)=TG 
DO  26  1=1 .NE 
Hljn)=0 
HJKn)=0 
HKI(I)-0 
TIJn)-0 
TJK(I)=0 
TKI(I)=0 
HHIJI 
HHJK( 

HHKII 
CONTINUE 

READING  IN  THE  ELEMENT  NUMBER  AND  ITS  NODAL  POINTS  AND  THE 
MATERIAL  TYPE.  WHICH  INCLUDES 

MAT(J)  = 1 SOIL  IN  INNER  EARTH  REGION 
= 2 PIPE  INSULATION 

=3  AIR  SPACE  SURROUNDING  THE  PIPE  IN  CASING 
= 4 OUTER  CASING  OF  THE  INSULATED  PIPE 
- 5 POURED-IN  INSULATION  MATERIAL  IN  BACK-FILL 
REGION 

= 6 NATIVE  SOIL  IN  THE  OUTER  EARTH  REGION 


DO  30  1=1 .NE 

READ(9 . •)  J . (NODE(J ,K) .K=1 .3) ,MAT(J) 


IF  (MAT( 

;j; 

i.EQ.i; 

1 CC(J 

l-KTCT/12 

IF  fMAT< 

J 

I .EQ.2 

1 CCfJ 

l-KI/12. 

IF  (MAT( 

J 

I.E0.3 

1 CC(J 

l-KASP/12 

IF  (MAT( 

J 

1 .E0.4 

1 CC(J 

l-KCAS/12 

IF  fMAT( 

J 

1 .EQ.5 

1 CC(J 

l-KBF/12, 

IF  (MAT( 

1 .EQ.6’ 

1 CC(J] 

l-KG/12. 

CONTINUE 

READ  IN  TOTAL  NUMBER  OF  ELEMENTS  HAVING  BOUNDARY  SEGMENTS  SUBJECT 
TO  CONVECTIVE  HEAT  TRANSFER 
READ  (9.*)  NECB 

C READ  IN  ELEMENT  NUMBER.  CONVECTIVE  HEAT  TRANSFER  COEFFICIENTS. 

C AND  AMBIENT  TEMPERATURES  FOR  THREE  BOUNDARY  SEGMENTS 
DO  35  1=1. NECB 

READ  (9.*)  M.HIJ(M),HJK(M).HKI(M).TIJ(M).TJK(M).TKI(M) 
IXCB(M)-1 


35 

CONTINUE 

ITER=1 

38 

DO 

40  1-1 .NN 

DO 

40  J=1 .NN 

Q(1 

[)=0. 

KK( 

;i.J)-0 

OQI 

I)-0. 

CK< 

I . J )-0 

DQ( 

!l)-0. 
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W(I)-=1  .0 
INDX(I)-:1 
40  CONTINUE 
L-1  . 

DO  180  1^1  .NE 
I-NODE(M, 1 ) 

J-=NOOE(M.2) 

K-NODE(M.3) 

IF(MAT(M) .EQ.3)  CC(M)-KASP/12 . 

C=CCfM) 

IF  ((INXK.EQ.0) .OR. (ITER.E0.1))  GO  TO  60 
C DETERMINE  SOIL  AND  INSULATION  THERMAL  CONDUCTIVITIES  BASED  ON  THE 
C MEAN  TEMPERATURES 

TM=(T(I)+T(J)+T(K))/3. 

IF(MATfM) .EQ.2)  CALL  INSULK(TM.C. ITYPE) 

IFfMATfM) .EQ.e)  CALL  SOILKfTM.KG.C) 
if(mat(m).eq.i)  call  SOILK(TM.KG.C) 

CC(M)=C 

60  xi«=xn) 

XJ=X(J) 

XK=X(k) 

yi=y(i) 

YJ=Y( J) 

YK«=Y(K) 

CXX-=C 
CXY-=0 . 

CYX«=0. 

CYY-C 

B2I-YJ-YK 

B3I-:XK-XJ 

B2J-YK-YI 

B3J-XI-XK 

B2K=YI-YJ 

B3K«XJ-XI 

C CALCULATE  THE  ELEMENT  AREA 

SA=0.5*(XJ*B2J+XI*B2I+XK*B2K) 

SA**ABS(SA) 

A2«=SA*2. 

AS(M)«A2 

B2I-B2I/A2 

B3I-B3I/A2 

B2J-=B2J/A2 

B3J=B3J/A2 

B2K-B2K/A2 

B3K-B3K/A2 

B2IZ(M)-B2I 

B3IZ(m)-B3I 

B2JZ(M)-B2J 

B3JZ(m)-B3J 

B2KZ(M)-B2K 

B3KZ(M)-B3K 

B I I-SA*  L* (B2 1 -BZ I •CXX+BZ I -BSI •CXY+B3I •BZ I •CYX+B3I •B3I •CYY) 

B I J-=SA*  L* ( B2 1 •BZ J *CXX+B2 1 •B3 J •CXY+B3 1 •BZ J •CYX+B3 1 -BSJ -CVY ) 

B I K-:SA*  L* ( BZ I -BZK^CXX+BZ I •B3K*CXY+B3 1 •B2K*CYX+B3 1 •BSK^CYY ) 

B J I-=SA*  L*  ( BZ  J •BZ  I •CXX+BZ  J •B3 1 •CXY+B3  J •BZ  I •CYX+B3  J •BS I •CYY) 
BJJ-=SA^L^(B2J^B2J^CXX+B2J^B3J^CXY+B3J^B2J^CYX+B3J^B3J^CYY) 
BJK-=SA^L^(B2J^B2K^CXX+B2J^B3K^CXY+B3J^B2K^CYX+B3J^B3K*CYY) 

BK I-SA^  L^ ( BZK^BZI •CXX+BZK^BSI •CXY+B3K^B2I •CYX+B3K^B3I •CYY) 
BKJ-:SA^L^IB2K^B2J^CXX+B2K^B3J^CXY+B3K^B2J^CYX+B3K^B3J^CYY) 
BKK-=SA^L^(B2K^B2K^CXX+B2K^B3K^CXY+B3K^B2K^CYX+B3K^B3K^CYY) 


KK( 

;i  .1] 

l-KK( 

;i . n 

l+BII 

KK( 

I .J 

l-KK( 

I .J 

l+BIJ 

KK< 

I .K 

l-KK( 

1 .K 

l+BIK 

KK( 

J.I 

l-KK( 

J.I 

l+BJI 

KK( 

J.J 

)-KK( 

J.J 

l+BJJ 

KKI 

J.K 

I -KKI 

J.K 

l+BJK 

KK( 

K.I 

l-KKI 

K.I 

f+BKI 

KKI 

K,  J 

>-KK| 

K.J 

l+BKJ 

KKI 

K.K 

I -KKI 

K.K 

l+BKK 

IFI 

[IXCB(M) , 

.EQ.0)  GO 

C ADDITION  OF  CONVECTION  TERMS  TO  THE  ELEMENT  MATRIX  TO  ACCOUNT 
C FOR  CONVECTION  ON  BOUNDARY 

C READING  IN  CONVECTIVE  HEAT  TRANSFER  COEFFICIENTS  AND  AMBIENT 
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C TEMPERATURES  FOR  THREE  BOUNDARY  SEGMENTS 

HHIJ(M)-HIJ(M)*L*SORTr fX(n-X(J))**2+(Y(I)-Y( J))**2)/6. 
HHJK(M)-HJKfM)*L*SQRT( (Xf j)-XfK))**2+(Y( j)-Y(K))**2)/6. 
HHKI(m)-HKI(M)*L*S0RT((X(K)-X(I))**2+(Y(K)-Y(I))**2)/6. 

KK( I . l)-=HHIJ(M}*2.+HHKI(M)*2.+KK(I .1) 

KKfl  . J )=HHIj(M)+KKn  , J) 

KKfl  .K)=HHKHM)+KKn  .K) 

KKf  J.  I)«=HHIJ(m)-H<K(J,  I) 

KKf J.J)=HHIJ(m)*2+HHJK(M)*2.+KK(J . J) 

KKf J ,k)=HHJK(m)+KK(J .K) 

KK(K, I )=HHKI (M)+KKfK, l) 

KKfK, J)-HHJK(m)+KK(K, J) 

KK(K.K)=HHJK(M)*2.+HHKI(M)*2.+KK(K.K) 

HHIJ(M)*=TIJfM)*3.*HHIJfM) 

HHJK(M)=TJKfM)*3.*HHJK(M) 

HHKI  (M)=TKI  (M) *3 . •HHKHm) 

130  OfI)-Q(n+HHIJfM)+HHKHM) 

Q(J)=Q(J)+HHIJ(m)+HHJK(m) 

Q(K)-0(KHHHJK(M)+HHKI  (M) 

180  CONTINUE 

185  FORMATf'  M I J K MAT.  C’) 

187  FORMAT(1X.5I6.F10.4) 

C DETERMINE  OUTER  SURFACE  TEMPERATURES  OF  UNDERGROUND  PIPES 
DO  200  1=1 .8 
T(I)=TP1 
11=1+8 
T(II)=TP2 
200  CONTINUE 

C DETERMINE  OUTER  BOUNDARY  TEMPERATURES  OF  EARTH  REGION 
CALL  TGXX(T.TGX. MONTH) 

MZ1=MZ+1 
DO  260  I=MZ1 .NN 
SUM=0. 

DO  250  J=1 .MZ 

250  SUM=SUM+KK(I.J)*T(J) 

QQ(I)-Q(I)-SUM 
260  CONTINUE 

IF(ICALB.EQ.I)  THEN 
WRITE(7.280) 

280  FORMAT (6X. *00  ARRAY') 

WRITE(7.285)  (QQ( I ) . 1=1 .NN) 

285  FORMAT  (5E12.5) 

END  IF 

C RENAMING  OF  MATRICES 
MNMNN-MZ 
DO  300  1=1 .MN 
K=MZ+I 

DO  290  J=1 ,MN 
KL-MZ+J 

290  CKf I . J)=KK(K.KL) 

XT(I)=T(K) 

DQ(I)=00(K) 

300  CONTINUE 

C SOLUTION  OF  SIMULTANEOUS  EQUATIONS 
C SET  PHYSICAL  DIMENSION  OF  MATRIX  A 
CALL  SOLVLE(CK.MN,NP,INDX,W,DQ) 

DO  310  1=1 .MN 
K-MZ+I 
T(K)«DQ(I) 

310  CONTINUE 

C CALCULATIONS  OF  EQUIVALENT  THERMAL  CONDUCTIVITIES  FOR  AIR 
C SPACE  INSIDE  THE  OUTER  CASING 
C 

C CALCULATE  AVERAGE  OUTER  SURFACE  TEMPERATURE  OF  INSULATED  PIPES 
SUl-0.0 
SU2=0 . 0 
DO  312  1-1 .8 
L1-I+NINS1-1 
L2-I+NINS2-1 
SU1-SU1+T(L1) 

SU2-SU2+T(L2) 

312  CONTINUE 

TSM1-SU1/8. 
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TSM2-SU2/8 . 

C DETERMINE  THE  EFFECTIVE  SURFACE  TEMPERATURE  OF  INSULATED  PIPES, 

C IN  DEG  F. 

DIP1-D1P+2.-THK1 

DIP2-D2P+2.*THK2 

TEFPS»(DIP1*TSM1+DIP2*TSM2)/(DIP1+DIP2) 

C CALCULATE  THE  INNER  SURFACE  TEMPERATURE  OF  THE  CASING 
SU3-0 . 0 
DO  314  1-1 ,8 
L3-I+NSIC-1 
314  SU3-SU3+T(L3) 

TSIOSU3/8. 

C DETERMINE  THE  TEMPERATURE  DIFFERENCE  BETWEEN  THE  EFFECTIVE  SURFACE 
C TEMPERATURE  OF  INSULATED  PIPES  AND  THE  INNER  SURFACE  OF  THE  CASING. 
C IN  DEG  F 

TDEL-ABS(TEFPS-TSIC) 

TAS=(TEFPS+TSIC)/2. 

CALL  EOUIKO(TAS,TDEL.EINS.ECAS.TEPS.TSIC.KASP) 

WRITE(7.320)  KASP 

320  FORMATf*  KASP-’ . FI 0 . 4 .2X . ’ (BTU-IN./H-FT**2-DEG  F) * ) 

330  FORMAT r TEMPERATURE  ARRAY  : T(I).  1=1 .NN  *) 

C CALCULATE  THE  MEAN  VALUES  OF  INSULATION  THERMAL  CONDUCTIVITY  FOR 
C PIPES  1 AND  2 
350  SK1 1-0. 

SKI 2-0. 

DO  400  LN-1 ,16 
LM-MINS+LN-1 
LL-LM+16 

SKI  1-SKI 1+CC(LM) 

SKI2-SKI2+CC(LL) 

400  CONTINUE 

KIX1-SKI1/16. 

KIX2-SKI2/16. 

R1-D1/24. 

R2-D2/24. 

TH1X-TH1/12. 

TH2X-TH2/12. 

IF(ICALB  .EQ.  0)  MO-11 

C CALCULATIONS  OF  THE  HEAT  LOSSES  FROM  THE  UNDERGROUND  PIPES 
CALL  PIPEHL(T.R1 .R2.TH1X.TH2X,KIX1 ,KIX2 ,MO,QTX) 

HLOSS-QTX 

IF(ITER.EO.I)  HLOSX-0. 

C DETERMINE  IF  PIPE  HEAT  LOSS  VALUE  HAS  CONVERGED,  OR  CONTINUE 
C ITERATIONS  IF  REQUIRED 

D E LQT-ABS ( H LOSS-H  LOSX ) /H  LOSS 

IF(DELOT.LE.  0.010)  GO  TO  2010 

ITER-ITER+1 

HLOSX-HLOSS 

GO  TO  38 

2000  WRITE  (7,2005) 

2005  FORMAT  (IX. 'THERE  ARE  SOME  ERRORS  IN  INPUT  DATA') 

2010  WRITE  (7.185) 

DO  2020  1-1 .NE 

WRITE  (7,187)  I.(NODE(I.J),J-1.3).MAT(I).CC(I) 

2020  CONTINUE 

WRITE  (7.330) 

WRITE  (7.285)  (T( I ) . 1-1 .NN) 

CALL  PIPEHL(T.R1 ,R2,TH1X.TH2X.KIX1 ,KIX2,7,QTX) 

STOP 

END 

SUBROUTINE  TGO(TGX.PI .Y) 

C THIS  SUBROUTINE  CALCULATES  THE  UNDISTURBED  EARTH  TEMPERATURES 
C AT  VARIOUS  DEPTHS 

DIMENSION  TGX(12.7) .Y(150) 

C READING  IN  THE  ANUAL  AVERAGE  TEMPERATURE  AND  AMPLITUDE  OF  THE 
C MONTHLY  NORMAL  TEMPERATURE  CYCLE  OF  THE  SITE,  IN  DEG  F,  AND 
C THERMAL  DIFFUSIVITY  OF  SOIL.  IN  FT**2/H. 

READ  (8.*)  AO.BO.DIFF 
W-2. •PI/12. 

WZ-2 . •PI/(8760*DI FF*2) 

Z2-SORT(WZ) 

DO  1 1-1,12 
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DO  1 J-1  ,7 
2-2Z*Y(33-J) 

TGX(I , J)-AO+BO*EXP(-2)*SIN(W*(I-3)-2) 

RETURN 
END 

SUBROUTINE  TGXX(T . TGX .MONTH) 

C THIS  SUBROUTINE  PROVIDES  OUTER  BOUNDARY  TEMPERATURES  OF  EARTH  REGION 
DIMENSION  T(150).TGX(12.7) 

T(32)=TGX(MONTH.1) 

DO  1 1=1 .8 
11=1+32 

1 T(II)=T(32) 

DO  5 1=2.7 
115=1+15 
J 1=33-1 

T(I15)=TGX(MONTH.I) 

T(JI)=TGX(MONTH.I) 

5 CONTINUE 

DO  10  1=1,3 
122=1+22 

10  T(I22)=TGX(MONTH.7) 

RETURN 

END 

SUBROUTINE  INSULK(TM.C. ITYPE) 

C THIS  SUBROUTINE  DETERMINES  THE  THERMAL  CONDUCTIVITY  OF  PIPE 
C INSULATION  (CALCIUM  SILICATE  AND  MINERAL  WOOL)  AS  A FUNCTION 
C OF  THE  MEAN  TEMPERATURE. 

REAL  KN(16).KINS 
DIMENSION  TN(16) 

DATA  KN  /0. 375. 0.40. 0.42. 0.45, 0.48. 0.50, 0.53, 0.555. 0.58.0. 61 . 

& 0.63.0.66.0.68.0.74,0.82.0.90/ 

C DETERMINE  THE  THERMAL  CONDUCTIVITY  OF  CALCIUM  SILICATE 
IF  (ITYPE  ,EQ,  1)  THEN 
DO  5 J=1 .16 
IF(J  .LE.  13)  THEN 

TN(J)=100.+(J-1)*50. 

ELSE 

TN(J)=700.+(J-13)*100. 

END  IF 

5 CONTINUE 

IF(TM  .GT.  TN(1))  GO  TO  10 
KINS-KN(I) 

GO  TO  100 

10  IF(TM  .LT.  TN(16))  GO  TO  20 

KINS-KN06) 

GO  TO  100 
20  DO  50  1-1,15 
TI-TM-TN(I) 

IF(T1  .NE.  0.)  GO  TO  30 
KINS-KN(I) 

GO  TO  100 

30  T2-TN(I+1)-TM 

IF(T2  .NE.  0.)  GO  TO  40 
KINS-KN(I+1) 

GO  TO  100 
A0  P-T1*T2 

IF(P  .LT.  0.)  GO  TO  50 

KINS-KN(I)+T1*(KN(I+1)-KN(I))/(TN(I+1)-TN(I)) 

GO  TO  100 
50  CONTINUE 
100  C-KINS/12. 

C DETERMINE  THE  THERMAL  CONDUCTIVITY  OF  MINERAL  WOOL 
ELSE  IF  (ITYPE  .EQ.  2)  THEN 

K INS-0. 2420  + TM*1.5501E  -4  + TM*TM*7.5001E  -7 
C-KINS/12. 

ELSE 
RETURN 
END  IF 
RETURN 
END 
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SUBROUTINE  SOI LK(TM. KG , C) 

C THIS  ROUTINE  DETERMINES  THE  THERMAL  CONDUCTIVITY  OF  SOIL  AS  A 
C FUNCTION  OF  MEAN  TEMPERATURES. 

REAL  K(14) .KG 
DIMENSION  TX(14) 

DATA  K/1 . 1 , 1 . 1 , 1 . 1 , 1 .0,0.4,0.31 .0.25,0. 19,0. 15,0. 1 1 ,0.09,0.07, 
ft  0.05,0.05/ 

DO  1 I-=1 .14 

1 TX(I)-50.+(I-1)*25. 

IF(TM.GT.TX(1 ))  GO  TO  5 
ZK=1 . 1 
GO  TO  50 

5 IF(TM. LT.TX(14))  GO  TO  10 

ZK-0.05 
GO  TO  50 

10  DO  25  I-l .13 
T1=TM-TX(n 
IF(T1 .NE.0)  GO  TO  15 
ZK=K(I) 

GO  TO  50 
15  CONTINUE 

T2=TM-TX(I+1) 

IF(T2.NE.0.)  GO  TO  20 
ZK=K(I+1) 

GO  TO  50 
20  CONTINUE 

P-:T1*T2 

IF(P.GT.0)  GO  TO  25 

ZK-K( 1+1 )+T2* (K( 1+1 )-K( I ) )/25 . 

GO  TO  50 
25  CONTINUE 

50  C-=ZK.KG/(1  . 1*12.  ) 

RETURN 

END 

SUBROUTINE  PIPEHL(T.R1 .R2.TH1 .TH2.ZKS1 ,ZKS2,MO.QT) 

C THIS  SUBROUTINE  CALCULATES  THE  AVERAGE  TEMPERATURE  DROPS  ACROSS  THE 
C PIPE  INSULATIONS  AND  THE  RATES  OF  HEAT  LOSS  FROM  THE  UNDERGROUND 
C PIPES  IN  DIRECTLY  BURIED  CONDUIT  SYSTEM 
DIMENSION  T(150) 

PI-4.*ATAN(1  .) 

SUM 1=0. 

SUM2=0 . 

N1=8 

DO  1 1=1 ,N1 
K1  = I 
K2=I+8 
K3-I+84 
K4=I+92 

SUM1=SUM1+T(K1 )-T(K3) 

SUM2=SUM2+T ( K2 ) -T ( K4 ) 

1 CONTINUE 

T1=SUM1/N1 
T2=SUM2/N1 
2KIS1=ZKS1*12. 

ZKIS2-ZKS2*12. 

Q1=ZKS1*2.*PI*T1/LOG((R1+TH1)/R1) 

Q2=ZKS2*2 . -PI •T2/LOG( (R2+TH2)/R2) 

QT=Q1+Q2 

IF(MO  .EO.  11)  GO  TO  50 
WRITE(M0.5)  ZKIS1.ZKIS2 

5 FORMAT(/’  AVERAGE  VALUES  OF  PIPE  INSULATION  THERMAL’ . 

ft’  CONDUCTIVITY  KI1  - ’,F10.3.’  KI2  - ’.F10.3, 

ft  ’ BTU-IN/H-FT**2-DEG  F ’) 

WRITE(MO.10)  T1 .T2 

10  FORMAT(/’  AVERAGE  TEMPERATURE  DROPS  ACROSS  INSULATION  : ’ ,/, 
ft  ’ T1«  •.F10.2,’  T2=  ’,F10.2.’  DEG  F’) 

WRITE(MO.20)  01 ,Q2,QT 

20  FORMAT (/.2X. ’HEAT  LOSSES  FROM  UNDERGROUND  PIPES  : ’/’  01=’. 

ft  F10.2.  ’ O2=’.F10.2.’  OT=’,F10.2.’  BTU/H-FT’) 

50  RETURN 

END 
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SUBROUT INE  PIPEJ(X.Y. TRTK . D , F . I NXK . I TYPE ) 

C THIS  SUBROUTINE  READS  IN  THE  INPUT  DATA  TO  BE  USED  FOR  CALCULATIONS 
C OF  THE  HEAT  LOSSES  FROM  THE  UNDERGROUND  PIPES  AND  GENERATES  X AND  Y- 
C COORDINATES  OF  NODAL  POINTS  FOR  THE  TWO  PIPE  SYSTEM. 

REAL  KII .KIG.KI .KG 
DIMENSION  X(150) .Y(150) 

COMvKDN  /PP/TP1 .TP2.KII .KIG.DII .DI2,THI1 .THI2.B1 ,B2.S1 .S2.TG. 
ft  WW.HY, MONTH 

COMADN  /EKC/D 1 P . D2P . D I C . THK 1 . THK2 
C READ  TEMPERATURE  OF  PIPE  NUMBERS  1 AND  2.  IN  DEG  F 
READ  (8,*)  TP1 .TP2 

READ  THERMAL  CONDUCTIVITY  OF  THERMAL  INSULATION  AND  SOIL. 
RESPECTIVELY.  IN  BTU-IN./H-FT.*2  - DEG  F.  AND  INDEX  OF  THERMAL 
CONDUCTIVITY  : INXK  = 0 CONSTANT  THERMAL  CONDUCTIVITY 

= 1 TEMPERATURE  DEPENDENT  THERMAL  CONDUCTIVITY 
AND  THE  MATERIAL  TYPE  OF  PIPE  INSULATION.  WHICH  INCLUDES 
ITYPE  «=  1 CALCIUM  SILICATE 
= 2 MINERAL  WOOL 
READ  (8.*)  KII .KIG. INXK. ITYPE 
C READING  IN  THE  OUTSIDE  DIAMETERS  OF  STEEL  PIPES  1 AND  2.  IN  INCHES 
READ  (8.*)  Dll .DI2 

C READING  IN  THE  THICKNESS  OF  THERMAL  INSULATION  USED  FOR  PIPES  1 
C AND  2.  RESPECTIVELY.  IN  INCHES 
READ  (8.*)  THI1 .THI2 

C READ  IN  THE  INSIDE  DIAMETER  AND  THICKNESS  OF  OUTER  CASING.  AND  THE 
C DEPTH  OF  ITS  CENTER  BELOW  GROUND  SURFACE.  IN  INCHES. 

READ  (8.*)  DIAC.THKC.DEPC 

C READING  IN  THE  VERTICAL  DISTANCES  (IN  FT.)  FROM  HORIZONTAL 

C CENTERLINE  OF  THE  INNER  EARTH  REGION  TO  CENTERS  OF  PIPES  1 AND 

C 2.  RESPECTIVELY.  AND  THE  AVERAGE  EARTH  TEMPERATURE.  IN  DEG  F. 

READ  (8.*)  SI .S2.TG 

C READ  IN  THE  THICKNESS  OF  EARTH  COVER. WIDTH  AND  DEPTH  OF  OUTER 

C EARTH  REGION  SURROUNDING  THE  UNDERGROUND  SYSTEM.  IN  FT. 


READ  (8.*)  E. WW.HY 

WRITE(7.10)  TP1 .TP2. KII .KIG.DII .DI2 
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FORMATC 

TP1  TP2  KI 

KG 

D1  D2*/6F7.2) 

WRITE(7.20) 

THI1 .THI2. DIAC.THKC.DEPC. SI .S2.TG 

20 

FORMAT ( ’ 
ft  /8F7.2) 

THI1  THI2  DIAC 

THKC 

DEPC  SI 

S2 

TG' 

WRITE(7.30) 

E. WW.HY. MONTH 

30 

FORMATC 

E WW  HY 

MONTH 

V3F7.2.I7) 

C READ  IN  THE  INSIDE  WIDTH  AND  HEIGHT  OF  THE  INNER  EARTH  REGION. 

C FOR  LOOSE-FILL  INSULATION  SYSTEM.  IN  FT. 

READ  (8.*)  A.B 

C CHANGE  TO  ENGINEERING  UNITS 
D1=DI1/12. 

R1-DU0.5 

D2-DI2/12. 

R2«=D2*0.5 

D1P»=DI1/12. 

D2P-DI2/12. 

DEP-DEPC/12. 

DIC-DIAC/12. 

KI-KII/12. 

KG=KIG/12. 

W=A+2*TRTK/12 

h^B+E+D+F 

WRITE(7.40)  W.H.D.F. A.B. WW.HY 

40  FORMATC  WH  D F A BWWHY’ 

ft./.8F7.2) 

PI«4.*ATAN(1.) 

TH1=THI1/12. 

TH2-THI2/12. 

THCA=THKC/12. 

THK1«=THI1/12. 

THK2=THI2/12. 

C DETERMINE  THE  X AND  Y-COORDINATES  OF  EARTH  COVER.  SIDES  AND  BOTTOM 

C BED  OF  THE  INNER  EARTH  REGION  (NODAL  POINTS  41  TO  72) 

DO  50  1-1. 5. 2 
I 67-1+67 

X(I67)-W-(I-1)*W/4. 

50  YU67)-E 

DO  60  1-1. 3. 2 
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I68-I+68 

X(I68V(W+A)*e.5  - A*(I-1)*0.5 
60  Y(I68)-E 

DO  65  1-1 .4 
140-1+40 
1 44—1+44 
148-1+48 
1 52- 1 +52 
I1-I-1 

X(I40}-fW-A)*0.5 

Y(I40)-(D+E)+I1*B/4. 

X(I44)-0.5*(W-A)+I1*A/4. 

Yf I44)-D+E+B 

X( I48)-fW+A)*0.5 

Y( 148)-(D+E+B)-I 1 *B/A. 

X( I52)-0.5*(W+A)-I1*A/4. 

YU52)-D+E 
65  CONTINUE 

DO  70  1-1 .3 
156=1+56 
160=1+60 
164=1+64 
X(I56)-0.0 

Y(l56)=(D+E)+B*(I-1)/2. 

Xf I60)=(W-A)*0.5+(I-1)*A*0.5 

Y(I60)-H 

Xn64)-W 

Y(I64)-(D+E)+(3-I)*B/2. 

70  CONTINUE 

Xf60)-0.0 

Yf60)-H 

Xf64)-W 

Y(64)-H 

C THE  X AND  Y-COORDINATES  OF  OUTER  BOUNDARY  EARTH  SURROUNDING  THE 
C DIRECTLY  BURIED  CONDUITS  (NODAL  POINTS  17  TO  40. AND  73  TO  84) 
X(17)— WW 
YM7)-E+0.083 
X(31 )-W+WW 
Y(3l)-E+0.083 
DO  72  1-1,2 
I17-I+17 
128=1+28 
Xf  117)— WW 

Yni7)-(O+E)+I*B*0.25 

X(I28)-W+WW 

Y(I28)-(D+E)+(3-I)*B*0.25 
72  CONTINUE 

DO  75  1-1,3 
I19-I+19 
I22-I+22 
125=1+25 
X(I19)— WW 

Yni9)-H+HY*(I-1)*0.5 
Xn22)-W.(I-1)*0.5 
Y(  I22)-h++HY 
X(I25)-W+WW 
Y(I25)-H+HY*(3-I)*0.5 
75  CONTINUE 

DO  77  1-1 .5.2 


I33-I+33 

X( 

;i33^ 

l-W-(I-1)*W/4. 

77  Y( 

133 

1-0.0 

DO  78  1-1 

1 .3.2 

1 34- 1 +34 

X( 

;i34; 

|-(W+A)*0.5-A*(I- 

78  Y( 

134 

1-0.0 

DO  80  1-1 

1 .2 

I31-I+31 

I38-I+38 

X( 

ri3i ; 

)-W+WW*(3-I)*0.5 

Yl 

131 

)— 0 . 

XI 

138 

I— WW*0.5*I 

Yl 

(138; 

1-0. 
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80 


82 


85 

C 


90 

C 

C 

C 


95 

C 

C 
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C 

C 

C 


C 

C 

C 


CONTINUE 
DO  82  1-1  .2 
I 72- 1+72 
174-1+74 
1 80- 1 +80 
1 82- 1 +82 
X(I72)— WW.0.5 
Yn72)-E+(I-1)*D 
X(I74)— WW*0.5 

Yn74)=(E+D+B*0.5)+(I-l)*(F+B*0.5) 

X(I80)-W+0.5*WW 

Yn80)  = (E+D+B+F)-(I-1)»(F+0.5*B) 

Xn82)-W+0.5*WW 
Y(I82)-(E+D)-(I-1)*D 
CONTINUE 
DO  85  1=1 .4 
176=1+76 

X(I76)=W*(I-1)/3. 

Y(I76)=H+HY*0.5 

X AND  Y-COORDINATES  OF  THE  CENTERS  OF  THE  PIPES 
XC1-W*0.5 
B1-DEP-S1 
YC1-B1 
XC2=W*0.5 
B2-DEP+S2 
YC2-B2 

WRITE(7.90)  XC1 .YC1 .XC2.YC2,B1 .B2 

FORMATC  XC1  YC1  XC2  YC2  B1  B2’/6F7.3) 
THE  X AND  Y-COORDINATES  OF  NODAL  POINTS  AT  THE  INNER  AND  OUTER 
SURFACES  OF  PIPE  INSULATION  (NODAL  POINTS  1 TO  16.  AND  85  TO 
100) 

DO  95  1=1 ,8 

THETA-2. •PI *1/8. 

18=1+8 
I 84- I +84 
192=1+92 

X( I )-XC1+0 .5*D1 •SIN(THETA) 

Y( I )-YC1+0 . 5*D1 •COS(THETA) 
xn8)=XC2+0.5*D2*SIN(THETA) 

Y f 1 8 ) -YC2+0 . 5*  D2 • COS ( THET A ) 


X( 

184' 

I-XC1+I 

;thi+ri ; 

••SIN! 

'THETA 

Y( 

184 

l=YC1+( 

TH1+R1 

••COS! 

THETA 

XI 

192 

I-XC2+I 

TH2+R2 

••SINI 

THETA 

Y( 

192 

l=YC2+( 

;th2+r2; 

••COSI 

[THETA 

CONTINUE 

THE  X AND  Y-COORDINATES  OF  NODAL  POINTS  AT  THE  INNER  AND  OUTER 
SURFACES  OF  PIPE  CASINGS  (NODAL  POINTS  101  TO  116) 

DO  97  1-1 .8 

THETA-2. •PI *1/8. 

1100=1+100 

1108=1+108 

Xri100)=XC1+DIC^0.5^SIN(THETA) 

Yni00)=DEP+DIC^0.5^COS(THETA) 

Xni08)=XC1+(THCA+DIC^0.5)^SIN(THETA) 

Y ( 1 1 08 )=DEP+ ( THCA+D I C^0 . 5 ) *008 ( THETA) 

CONTINUE 

THE  X AND  Y-COORDINATES  OF  NODAL  POINTS  AT  THE  CENTER  OF 
CONDUIT  CASING. AND  AT  THE  TOP  AND  THE  BOTTOM  OF  AIRSPACE 
(NODAL  POINTS  117  TO  119). 

X(119)=W^0.5 

Yf119)=DEP 

X(117)=W^0.5 

Y(117)-0.5^(Y(88)-Y(104)) 

X(118)=W^0.5 

Y(1 18)=0.5^(Y(108)-Y(100)) 

THE  X AND  Y-COORDINATES  OF  NODAL  POINTS  IN  BACK-FILL  SOIL.  OR 
POURED-IN  INSULATION  SURROUNDING  THE  PIPES  (NODAL  POINTS  120  TO 
131) 

YUP-0.5^(Y(112)-Y(55)) 

XLT-0.5^(X(114)-X(43)) 

YLO=0.5^(Y(47)-Y(116)) 

XRT-0.5^(X(51)-X(110)) 

DO  100  1-1 ,3 
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I119-I+119 
1 122-1  + 122 
I125-I+125 
I128-I+128 

X(1119)-0.5*(W+A)-e.25*A*I 
Y(I119)-([>+E)+YUP 
Xni22)-0.5*(W-A)+XLT 
Yf I122)-(D+E)+0.25*B*I 
X(I125)-0.5*(W-A)+0.25*A*I 
Yni25)-=(D+E+B)-YL0 
X (1 1 28 )-0 . 5* (W+A)-XRT 
Y(I128)-(D+E+B)-0.25*B*I 
100  CONTINUE 
RETURN 
END 

SUBROUTINE  TWDPIP( IREPT) 

C THIS  SUBROUTINE  DETERMINES  THE  HEAT  LOSSES  FROM  TWO  PIPES  TO  THE 
C UNDERGROUND  SURROUNDING  THE  HEAT  DISTRIBUTION  SYSTEM. 

REAL  KII .KIG 

COMMON  /PP/T1 .T2,KII .KIG.DII .DI2.THI1 .THI2.B1 .B2.S1 .S2,TG. 
k WW,HY, MONTH 
PI=4..ATAN(1 .) 

X1=2. ‘PI 
R1-DI1/24. 

R2=DI2/24. 

TH1X-THI1/12. 

TH2-THI2/12. 

ZK1-KII/12. 

ZK2-ZK1 

D1-B1 

D2-B2 

ZKS-KIG/12. 

DO  10  1-1 . IREPT 
2 TH1=TH1X+0.1*(I-1) 

S-S1+S2 

A-R 1 +R2+TH 1 +TH2+0 . 05 
THI1=TH1*12. 

IF(A  -LT.  S)  A-S 
C1=X1*ZK1/LOG((R1+TH1)/R1) 

C2=X 1 • ZK2/ LOG ( ( R2+TH2 ) /R2 ) 

P1 1-1 .+C1/(X1*ZKS)*LOG((2*D1)/(R1+TH1)) 

P12-C2/(X1*ZKS)*LOG((A*A+(D1+D2)**2)/(A*A+(D1-D2)**2))*0.5 
P21-C1/(X1*ZKS)*LOG((A*A+(D1+D2)**2)/(A*A+(D1-D2)**2))*0.5 
P22-1 .+C2/(X1*ZKS)*LOG((2*D2)/(R2+TH2)) 

DEL-P12*P21-P11*P22 
ZKP1-C1*(P12-P22)/DEL 
ZKP2-C2*(P21-P11)/DEL 
TP1-(P12*T2-P22*T1 )/(P12-P22) 

TP2-(P21 •TI-PI 1 •T2)/(P21-P1 1 ) 

01-ZKP1*(TP1-TG) 

Q2=ZKP2*(TP2-TG) 

OT-Q1+02 

TAVG-(T1+T2)*0.5 

ZK=OT/(TAVG-TG) 

WRITE(7.6)  Dll .DI2.S1 ,S2.THI 1 ,KI I .KIG.T1 .T2 
6 FORMATC  Dll  DI2  SI  S2  THI1  KII  KIG  TP1  TP2' . 
ft/.7F6.2. 1X,2F6.0) 

WRITE(7.8)  Q1.Q2.QT.ZK 

8 FORMATC  Q1  Q2  QT  KP*/.3F7.2.2X.F6.3/) 

10  CONTINUE 

RETURN 
END 

SUBROUT I NE  EOU I KO ( T AS . TD E L . E I NS . ECAS . T EPS . TS I C . KASP ) 

C THIS  ROUTINE  CALCULATES  EQUIVALENT  THERMAL  CONDUCTIVITY  OF  AIR 
C SPACE  SURROUNDING  INSULATED  PIPES  INSIDE  THE  OUTER  CASING. 

REAL  KASP 

C0M40N  /EKC/D1P.D2P,DIC.THK1  ,THK2 
PI-4.*ATAN(1 .) 

C CALCULATE  THERMAL  CONDUCTIVITY.  IN  BTU-FT/H-FT**2-DEG  F.  AND 
C KINEMATIC  VISCOSITY.  IN  FT**2/S.  OF  AIR 
THKAIR-0. 01319  + TAS*2.5E  -5 
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VAIR-1.2624E  -4  + TAS*5.4E  -7 

C CALCULATE  THE  EFFECTIVE  DIAMETER  OF  THE  INSULATED  PIPES. 

C IN  FT.  AND  THE  CHARACTERISTIC  LENGTH  OF  AIR  SPACE 
DEFPIP-(D1P+2 . -THKI )+(D2P+2. •THK2) 

CL-=(DIC  - DEFPIP)*0.5 

C CALCULATE  THE  PRANDTL  NUMBER  OF  AIR  . AND  GRASHOF  NUMBER.  AND 
C EQUIVALENT  THERMAL  CONDUCTIVITY  OF  AIRSPACE  DUE  TO  FREE 
C CONVECTION.  IN  BTU-IN/H-FT**2-DEG  F. 

C 

CLODEFPIP 

PR ANT L-0. 7 1849  - TAS  • 1.275E  -4 

GRASOF=32.2  • TDEL  •(CLC*»3. )/((VAIR**2. )*(TAS+459.7)) 

CONKA=12.  •THKAIR*0.530*  ( (PRANTL*GRASOF) **0 . 250) 

C CALCULATE  THE  EMISSIVITY  TERM  AND  EQUIVALENT  THERMAL 
C CONDUCTIVITY  OF  AIRSPACE  DUE  TO  RADIATIVE  TRANSFER.  IN 

C BTU-IN/H-FT**2-DEG  F. 

etot=i./e:ins  + (defpip/dic)*(i  ./ecas  --i.) 

TREP=(TEPS  + 459.7)/100. 

TRSI«=(TSIC  + 459.7)/100. 

HRAD=0. 1714/100.  • (TREP+TRSI)*(TREP*TREP  + TRSI •TRSI )/ETOT 
RADKA*12  0 • HRAD  •CL 

C CALCULATE  EQUIVALENT  THERMAL  CONDUCTIVITY  OF  AIRSPACE.  IN 

C BTU-IN/H-FT**2-DEG  F. 

KASP=CONKA  + RADKA 
WRITE(7.50)  PRANTL.GRASOF 
50  FORMATC  PRANTL=' .F10.4.2X. •GRASOF=’ .E12.4) 

WRITE(7.60)  CONKA. RADKA 

60  FORMATC  CONKA=’ .F10.4.2X. ’ RADKA=' .F10.4.2X. 

4’ (BTU-IN/H-FT**2-DEF  F)') 

RETURN 

END 

SUBROUTINE  SOLVLE(A.N.NP,  INDX.W.B) 

C GIVEN  AN  NXN  MATRIX  A.  WITH  PHYSICAL  DIMENSION  NP,  THIS  ROUTINE 
C REPLACE  IT  BY  THE  LU  DECOMPOSITION  OF  A ROWWISE  PERMUTATION  OF 
C ITSELF.  INDX  IS  AN  OUTPUT  VECTOR  WHICH  RECORD  THE  ROW  PERMUTATION 

C EFFECTED  BY  THE  PARTIAL  PIVOTING;  W IS  VECTOR  OF  SCALING  FACTORS. 

C 

C THIS  ROUTINE  IS  USED  TO  SOLVE  THE  LINEAR  SET  OF  EQUATIONS  : 

C [A][X]»[B] 

DIMENSION  A(NP.NP) . INDX(N) .VV(N) .B(N) 

FORM  IMPLICIT  SCALING  VECTOR  W 

DO  12  I«1.N 
AAMAX  ■=  0.0 
DO  11  J-1.N 

IF(ABS(A(I . J)).GT.AAMAX)  AAMAX»ABS(A( I , J ) ) 

11  CONTINUE 
IF(AAMAX.EO.0.)  THEN 

WRITE(7.100)  I 

00  FORMAT( IX. 'ERROR: SINGULAR  MATRIX  - ZERO  ROW  : ROW', 15) 

RETURN 
END  IF 

VV(I)  - 1.0/AAMAX 

12  CONTINUE 

CROUT  METHOD:  LOOP  OVER  COLUMNS 

DO  19  J-:1  .N 
DO  14  1=1 . J-1 
SUM  - A(I ,J) 

DO  13  K=1 . 1-1 

SUM  - SUM  - A(I .K)*A(K.J) 

13  CONTINUE 
A(I.J)  = SUM 

14  CONTINUE 

PIVOT  IMPLEMENTATION 

AAMAX  = 0.0D0 
DO  16  I=J.N 
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SUM  - A(I ,J) 

DO  15  K-1 . J-1 

SUM  - SUM  - A(I .K)*A(K, J) 

15  CONTINUE 
A(I.J)  - SUM 

DUM  - W(I)*ABS(SUM) 

IF(DUM.GE.AAMAX)  THEN 
IMAX  > I 
AAMAX  - DUM 
END  IF 

16  CONTINUE 

IF(J.NE. IMAX)  THEN 
DO  17  K*1 ,N 

DUM  -=  A(IMAX,K) 

AfIMAX.K)  -=  A(J,K) 

A(J,K)  - DUM 

17  CONTINUE 
W(IMAX)  * W(J) 

END  IF 

INDX(J)  - IMAX 
IF(A(J.J).EO.0.0)  THEN 
WRITE(7,110)  J 

110  FORMAT nx. 'ERROR:  SINGULAR  MATRIX  - ZERO  " DIAG  " : ROWM5) 

RETURN 
END  IF 

IF(J.NE.N)  THEN 
DUM  « 1 .0/A(J.J) 

DO  18  I-J+1 .N 

A(I.J)  - A(I.J)*DUM 

18  CONTINUE 
END  IF 

19  CONTINUE 

FORWARD  SUBSTITUTION 

II  -=  0 
DO  22  I»1.N 
LL  - INDX(I) 

SUM  - B(LL) 

B(LL)  - B(I) 

IF(II.NE.0)  THEN 
DO  21  J-II.I-1 

SUM  - SUM  - A(I ,J)*B(J) 

21  CONTINUE 
ELSE  IF(SUM.NE.0.0)  THEN 

II  - I 
END  IF 
B(I)  - SUM 

22  CONTINUE 

BACKWARD  SUBSTITUTION 

DO  24  I-N.1 ,-1 
SUM  - B(n 
IF(I.LT.N)  THEN 
DO  23  J-I+1 .N 
SUM  - SUM  - A(I . J)*B(J) 

23  CONTINUE 
END  IF 

B(I)  - SUM/A(I.I) 

24  CONTINUE 
RETURN 
END 
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